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Sommario 



La Stabilita Idrodinamica e la transizione alia turbolenza sono stati oggetto di 
studio sin dalla fine del XIX secolo. In particolare, la discrepanza tra teoria e 
osservazioni sperimentali nel caso di transizioni subcritiche, ha costituito un pro- 
blenia complicato che ha promosso la ricerca di altri meccanismi che potessero 
generare la transizione, differenti da quello classico che prevede la crescita espo- 
nenziale asintotica delle onde di Tollmien-Schlichting. Nonostante il ruolo delle 
nonlinearita sia universalmente riconosciuto, un rinnovato interesse verso I'analisi 
lineare a partire dalla fine del XIX secolo e derivato dai risultati dell'analisi non 
niodale. La possibilita di una crescita di tipo algebrico, pur significativa e anche 
per perturbazioni asintoticamente stabili come nel caso del fiusso piano di Couette, 
ha aperto un nuovo scenario nello studio sulla transizione laminare-turbolento. Si 
osserva infatti che alcuni meccanismi, come il vortex tilting oppure il ruolo delle 
perturbazioni ortogonali al fiusso medio, vengono riscontrati gia dall'analisi lineare 
e tridimensionale. Inoltre si e mostrato che la crescita in energia cinetica della 
perturbazione e soltanto attribuibile al un meccanismo lineare. 

Scopo del presente lavoro e quello di contribuire alle conoscenze attuali sull'e- 
voluzione temporale di piccole perturbazioni tridimensinali in fiussi confinati, in 
particolare il fiusso di Couette piano, tramite I'analisi delle velocita di fase e delle 
frequenze. La loro evoluzione temporale e stata, ed e tuttora, poco analizzata ma 
contiene in realta preziose informazioni sulla vita delle perturbazioni. I risultati 
ottenuti per fiussi di Couette e Poiseuille mostrano la possibilita di velocita di fase 
diverse per le tre componenti di velocita, e soprattutto la presenza di brusche varia- 
zioni o salti nell 'evoluzione temporale delle frequenze. Tali variazioni permettono 
di distinguere tre periodi distinti della vita della singola onda, V Early transient, 
V Intermediate transient e il Far transient, e sembrano essere correlate con I'instau- 
rarsi di certe condizioni di self-similarita nei profili di velocita o vorticita. Tali 
analisi non sarebbero state possibili senza lo sviluppo di un codice di calcolo in 
ambiente Matlab® basato su una soluzione semi analitica (per fiussi confinati) del 
problema ai valori iniziali di Orr-Sommerfeld e Squire. Tale soluzione e espressa 
come serie di funzioni ortogonali, e la soluzione approssimata viene ottenuta appli- 
cando il metodo variazionale di Galerkin. II codice risultante risulta decisamente 
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vantaggioso i termini di tempi di calcolo e accuratezza. A concludere il lavoro, vie- 
ne mostrata I'evoluzione lineare di disturbi localizzati in forma di pacchetti d'onda 
per flusso di Couette e di Strato Limite. Si evidenziano le analogie con uno scenario 
di transizione in presenza di spot turbolenti, le quali portano a supporre che alcune 
proprieta di tali strutture risiedano gia nelle equazioni di governo linearizzate. 



II presente lavoro di tesi e stato in parte svolto al dipartimento di Matematica del 
Massachusetts Institute of Technology di Cambridge (USA), sotto la supervisione 
della Prof.ssa Gigliola Staffilani e tramite il progetto di mobilita extra-UE FP 
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Chapter 1 



Introduction 



1.1 Linear Stability and transition 

The reasons for the breakdown of a laminar flow to turbulence has been one of 
the central issues in fluid mechanics for over a hundred years, for the many ap- 
plications in the engineering, meteorology, oceanography and astrophysics. The 
theoretical work on transition is mainly based on the linear stability studies, which 
were firstly initiated in the nineteenth century by Helmholts (1868), Rayleigh and 
Kelvin. Reynolds (1883) dedicated to experiments on the instability of the pipe 
flow, and was the first to find the existence of a critical velocity (actually, the non- 
dimensional parameter that now brings his name, the Reynolds number) above 
which the transition to turbulence occurs. He observed the intermittent character 
of this phase as well, naming flashes the objects that we now call turbulent spots. 
The formulation for the viscous stability problem is due to Orr (1907) and Sommer- 
feld (1908), who dedicated respectively to the Plane Couette flow and to the Plane 
Poiseuille flow. The Orr-Sommerfeld equation has become the basis of the modal 
theory of hydrodynamic stability. Many years later Tollmien (1929) calculated 
the first neutral eigenvalues for Plane Poiseuille flow, and Schlichting continued 
his work, leading to the deflnition of the TS-waves, whose role in the transition 
process is salient. 

Only in the second half of the twentieth century the three-dimensional initial value 
problem was considered. The transient dynamics of perturbations revealed aspects 
that made the non-modal problem even more of interest than than the past analysis 
on the asymptotic states. The most important result is the presence of an alge- 
braic behavior in the early and intermediate stages of a perturbation's life; three 
main reasons for the transient growth were found: the non-orthogonality of the 
eigenfunctions, the possible resonance between the Orr-Sommerfeld and the Squire 
solutions and, for unbounded or semi-bounded flows, the presence of a continuous 
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spectrum (see e.g. the works by Criminale and Gustavsson). The role of these 
mechanisms, though hnear, in a transitional scenario is evident, and it is easy to 
understand why many efforts were made in the last two decades to investigate the 
conditions for "optimal growth". Only in the recent years the role of the linear 
mechanisms in the subcritical transition to turbulence has been pointed out by 
many authors (see, among others, Henningson). 



1.2 Thesis motivations and layout 

The aim of the present thesis is to contribute to the actual knowledge about the 
transient behavior of small perturbations in channel flows. The focus will be on a 
quantity whose temporal evolution had not been considered in detail before: the 
phase velocity or, equivalently, the frequency of the components of velocity and 
vorticity of a perturbation. Throughout the present work, it will be shown as from 
the analysis of the wave frequency, three terms of a disturbance's life can clearly 
be discerned. Some properties of similarity of the velocity and vorticity proflles 
will also be highlighted. 

In Chapter 2 the mathematical background is given, and the principal equations 
and definitions are introduced. In Chapter 3 an analytical method to solve the 
Orr-Sommerfeld and Squire initial value problem is presented, together with the 
implementation of a Matlab® code to obtain approximate solutions. The suggested 
method is verified and used for the further analysis. The focus of Chapter 4 is on 
the perturbation frequency and phase velocity. Numerical results are shown in 
terms of both the vorticity and velocity components, and similarity properties of 
the profiles are investigated. The last Chapter concerns the evolution of wave 
packets and linear spots. 



1.2 - Thesis motivations and layout 




Figure 1.1: Turbulent spot on a flat plate. Top: Re = 200000, the sublayer of the spot is shown 
from the glass wall of the water channel, by a suspension of aluminium flakes (Cantwell et ai, 
1978). Bottom: Re — 400000, the typical arrowshape angle becomes steeper; visualization by 
smoke in air with flood lighting, photograph by R. E. Falco, taken from the book of Dyke (1982) 



Chapter 2 



Mathematical background 



2.1 Initial value problem for shear flows: viscous linear 
analysis 

2.1.1 Base governing equations 

In the present analysis, the flow is taken to be incompressible and the governing 
equations for infinitesimal disturbances in parallel fiows are considered. The base 
fiow general expression is Ui = U{y)6ii, i.e. the streamwise direction is x, and it 
only depends on the wall-normal direction y. The origin of the reference system is 
set on the channel symmetry plane xz for Plane Couette fiow and Plane Poiseuille 
fiow (PCf and PPf, in the following), and on the wall for Blasius boundary layer 
fiow (Bbl), i.e. the fiow along a fiat plate with zero pressure gradient(Fig. 2.1). 
The equations governing the general evolution of fiuid fiow are the Navier-Stokes 
equations, that using Cartesian tensor notation read 



dui dui dp 1 2 
at oxj oxi Re 


(2.1) 


OXi 


(2.2) 


supported with the typical initial and boundary conditions of the form 




w,(x,0) = M°(x) 

,. /,. i\ n „„ 11„ 


(2.3) 



■Ui(x, t) = on walls 

The physical quantities u,v and w represent the velocity components, p represents 
the fiow static pressure, and they appear in the system (2.1) in nondimensional 
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form. For PCf and PPf the reference length is the channel semi-height h, the 
reference velocity is assumed to be the medium wall velocity Up = {U~^ — U~)/2 
for PCf, and the centeline velocity Ucl ioT PPf. For Bbl, the velocity scale is the 
freestream velocity Uoo and the length scale is the boundary layer displacement 
thickness 6*, which takes the following expression, as exact solution of the Blasius 
equation (Schhchting, 1979, p. 141): 



ux 



6* = 1.7208^— — {displacement thickness) 



(2.4) 



The approximate expression for the geometric thickness, defined as the distance 

(2.5) 



for which u = 0.99f/oo, is found to be 



ux 



So.99 = 4.91./—— {geometric thickness) 
So the following definitions for the Reynolds number will be considered 

Re = — ^ Plane Couette flow Re = Plane Poiseuille flow (2.6) 



u 



Re 



U^6* 



V 



Blasius boundary layer 



(2.7) 



where v is the kinematic viscosity. The evolution equation for the disturbances 
can be obtained by splitting the flow in two components, the Base flow {Ui{y), P{y)) 
and the perturbed state ({tj(x, t),p(x, t)) so that the complete fluid fleld can be 
written as Mj = f/j + Ui and p = P + p. The nonlinear disturbance equations read 



duj 
dt 
duj 
dxi 



^U — -u-— - — + —V^u- - u 
'' dxi ■' dxn dxj Re 



duj 
^dxj 



(2.8) 
(2.9) 



toghether with the appropriate initial and boundary conditions. 
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(a) PCf: Lref = K Uref = Up (b) PPf: Uef = h, Uref = UcL (c) Bbl: Lref = S* , Uref = U^ 

Figure 2.1: Flow schemes: reference quantities, base flow and reference systems. 



2.1 - Initial value problem for shear Rows: viscous linear analysis 

2.1.2 Linearized perturbative equations 

Considering the x reference axis oriented as the base flow streamwise direction, 
so that it assumes the general expression f/j = U{y)6ii, the complete velocity fleld 
becomes u = (f/ + u,v,w). In particular, for PCf U{y) = y, and for Bbl the base 
velocity proflle is tabulated in self-similar coordinates (Rosenhead, 1963, Chap. 
V). Introducing the mean velocity proflle and assuming small perturbations, the 
following linear equations can be written, as shown by Schniid & Henningson (2001) 
and Criminale (2003): 



du 
dx 


dv 


f =0 

dz 




du 

dt 


dx 


.dU 
dy 


dp 1 „9~ 
dx Re 


dv 
dt 


dx 


dp 
dy 


Re 


dw 
dt 


dx 


dp 

dz 


Re 



(2.10) 
(2.11) 
(2.12) 
(2.13) 



Taking the divergence of the linearized momentum equations (2.11), (2.12), (2.13), 
and using the continuity equation (2.10), an equation for the fluctuating pressure 
can be obtained and used to eliminate the pressure terms, in combination with 
(2.12), leading to the following equation for the wall- normal velocity: 



"dt dx dy'^ dx Re 







(2.14) 



To completely describe the three-dimensional flow fleld, a second equation is nec- 
essary, and it is convenient to write an equation for wall-normal vorticity, deflned 
as fj 



du 
dz 


dw 
dx 


' d 
dt 



u 



d 



V 



dUdv 
dy dz 



dx Re 
The quantity F is then deflned as F = V'^v, so that the system becomes 



(2.15) 



h U V 

dt dx Re 



d 



1 



— + U—-—V^ 
dt dx Re 



W^v 
d^Udv 



dy"^ dx 
^ __dUdv 
dy dz 



(2.16) 
(2.17) 

(2.18) 



The perturbations are Fourier transformed in x and z directions: two real wavenum- 
bers, a and /3 are introduced along the x and z coordinates, respectively. The 
generic quantity / is hence expressed as 



f{y,t;a,(3) 



+ 00 



+ 00 



f{x,y,z,t)e-'''''''f^'dxdz 



(2.19) 



2 - Mathematical background 



The system can now be written in the following form 

- *^6 = f (2.20) 

1 = -"'cosmr + ^tcos(4./p + ^ (g - *^f ) (2.21) 
— = -ikcos{(P)Uri - iksin{4>) -^v + — (^— - k^j (2.22) 



where = tan~^{(3/a) is the perturbation obliquity angle, and k = \/ a^ + fP is 
the polar wavenumber. The following boundary conditons applies in the wavenum- 
ber space, respectively for Bbl and PCf: 

dv 
v{y = ±l,t) = —{y = ±l,t) = fiiy = ±l,t) = (2.23) 

dv 

v{y -^ +oo,t) = —{y -> +oo,t) = r}(y ^ +oo,t) = (2.24) 

The streamwise velocity u and the spanvise velocity w can be recovered from the 
following expressions 

^=-(a--/3r)) (2.25) 

^ = -(/3- + ar}) (2.26) 

2.1.3 Energy amplification factor 

In order to quantify the growth of the perturbations, a natural choice is the 
kinetic energy density, defined as 

1 P^ 

e{t;a,/3) = - {\u\^ + \v\^ + \w\'^)dy (2.27) 

2 Jya 



I rvb 
2 



Va 



dv 



dy 



+ k^\v\^ + \fiY]dy (2.28) 



where ya and 1/5 are the limits of the domain. As a disturbance measure, the proper 
quantity is the energy amplification factor, G, defined as the kinetic energy density 
normalized with respect to its initial value (Criminale et ai, 1997; Lasseigne et ai, 
1999) 

the temporal growth rate of the kinetic energy r is then introduced to evaluate the 
beginning of the exponential asymptotic period, when dr/dt — )■ 

u /Q^ log\e{t;a,/3)\ 
r{t; a, /3) = — , t > (2.30) 



Chapter 3 



Wave transient analysis: an eigenfunction 
expansion solution method 



3.1 Introduction 

In the present chapter an analytical solution to the Orr-Sommerfeld and Squire 
initial value problem (eq. 3.1 and 3.2) is researched for channel flows, aiming to a 
better understanding of the early and intermediate terms of a perturbation's life. 
As a starting point, the IVP in the normal- velocity and normal- vorticity form is 
considered 

(3.1) 



v{y = ±l,t) = —{y = ±l,t) = fi{y = ±l,t) = (3.3) 

v{y, t = 0) = vo{y) v{y, t = 0) = f,^{y) (3.4) 

where the prime symbol indicates a total derivative along y. The evolution of 
the wall-normal velocity v is described by the Orr-Sommerfeld PDE (3.1), which 
is of fourth order in the spatial coordinate y and homogeneous, with homogeneous 
boundary conditions. The Squire equation (3.2) is inhomogeneous and the forcing 
term — ■i/3f/(?/)'-0 is known as vortex tilting^ being the product of the main vorticity 
in the spanwise direction [Viz = —U') and the perturbation velocity v. This term 
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3 - Wave transient analysis: an eigenfunction expansion solution method 

is responsible of the increase of the normal vorticity, for three-dimensional pertur- 
bations (see Criminale et ai, 1997). 
About the initial conditions, the following will be used in the present work 



{y) = {l-yy Uy)=0 


Symmetrical 


(3.5) 


{y)=y{l-yy fjo{y) = 


Anti symmetrical 


(3.6) 



It is known that for bounded flows all eigenvalues of the Orr-Sommerfeld and 
Squire ODE are discrete and inflnite in number and that the eigensolutions of 
the problem form a complete set as proved by Schensted (1960) and DiPrima 
& Habetler (1969). For unbounded or semi-bounded flows (as the Wake or the 
Boundary layer flows) Miklavcic & Williams (1982) and Miklavcic (1983) proved 
that if the base flow decays in an exponential way, then only a flnite number of 
eigenvalues exists and a continuum is present, while if the decay is algebraic there 
exists a inflnite discrete set (without the continuum). 

The focus of this chapter is on channel flows. Most of the studies in the past cen- 
tury deal with the modal analysis. About the Orr-Sommerfeld ODE, it is possible 
to express the solution as a generalized Fourier series once a base of orthogonal 
functions is found, and variational or Galerkin methods can be applied to provide 
very accurate approximations when a flnite number of trial functions are used. 
The Orr-Sommerfeld and Squire modes can be used to express the solution (see 
Schmid & Henningson, 2001), however it was shown that some sets of normal func- 
tions can give better results in terms of accuracy and computational cost. Orszag 
(1971) solved the Orr-Sommerfeld ODE numerically using expansions in Chebyshev 
polynomials and used the Lanczos's tau method to determine the series coefficients. 
He showed that this series gives the highest convergence rate, since the error after 
N terms is smaller than any power of A^~^. Before him, Dolph & Lewis (1958) 
were the first applying a Galerkin method to obtain the coefficients (reduction to 
a system of A^ algebraic equations), together with the QR algorithm. They used 
normal functions that guarantee a A^^ rate. Gallagher & Mercer (1962) used the 
Chandrasekhar functions, adopted in the present work as well, which provide a 
rate of convergence of A^^. 

About the solution of the initial value problem, there is no conceptual difficulty 
in using the eigensolutions of the Orr-Sommerfeld and Squire ODE system but, as 
outlined in Drazin & Reid (2004), this requires the solution of the adjoint differ- 
ential equation. 

In the flrst part of this chapter an eigenfunction expansion method for the initial 
value problem (3.1)-(3.2) is proposed; the method does not involve the eigensolu- 
tions to the Orr-Sommerfeld and Squire ODE system, and the approximate time- 
dependent coefficients are obtained with the variational minimization principle. 
The two PDEs are then reduced to a system of A^ ODEs. A Matlab® code is 
implemented and verified and afterwards used for the analysis of Chapter 4, whose 
focus will be on the wave frequency and velocity profiles. 
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3.2 - Solution to v equation 



3.2 Solution to v equation 



3.2.1 Choice of a base of orthogonal functions 

The solution of (3.1) can be expressed as a generalized Fourier expansion, with 
time- dependent coefFcients: 



v{y,t) = Y,Cn{t)Xn{y) ye [-1,1] 



(3.7) 



ra=l 



where Xn{y) are orthogonal functions, and the following inverse transform applies 
(see Strauss, 1992): 



c„(t) 



J^^v{y,t)Xn{y)dy 
/_\x„(t/)X„(|/)d|/ 



(3.^ 



Since in the initial value problem both the initial condition and the boundary con- 
ditions need to be imposed, it is worthwhile to consider functions that satisfy the 
boundary conditions of the problem considered. Moreover, note that the coeffi- 
cients c„ of the series are in general complex, since v is complex-valued and the 
spatial modes are considered as real. The particular orthogonal functions which we 
use are those defined by the following fourth order eigenvalue problem satisfying 
the same boundary conditions of the original equation. This choice for the simpli- 
fied problem is not the only possible but revealed to be appropriate; note that this 
model equation is contained in the diffusive part of the PDE (3.1) 



dy^ 
X{y = ±1] 



\'X{y) y E [-1,1] 
dX , 

-0 ^iy-±i)-o 



(3.9) 
(3.10) 



A solution to this problem is obtained considering sines, cosines, hyperbolic sines 
and hyperbolic cosines (see appendix A.l for the complete solution). Two different 
sets of eigenvalues and the corresponding eigenf unctions are found, respectively 
odd and even, by numerically solving the following transcendental equations 



tan{\n) — tanh{\n) = {odd set) 
tan{\n) + tanh{\n) = {even set) 

The corresponding normalized eigenfunctions (Fig. 3.1) are 

'sinh{\ny) sin{\ny) 



Xr, 



Xn 




1 

71 



a/2 [ cosh{\n 



n = 1,3,5.., A^ — 1 {odd set) 
n = 2,4,6.., A^ {even set) 



(3.11) 
(3.12) 



(3.13) 
(3.14) 
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3 - Wave transient analysis: an eigenfunction expansion solution method 




Figure 3.1: The basis eigenfunctions 

Similar functions, in a different domain, have been used by Chandrasekhar (1961, 
app. V), in the study of the circular Couette flow between coaxial cylinders, and 
by Gallagher & Mercer (1962) to solve the Orr-Sommerfeld ODE. 

Since the imaginary and the real part of the solution -0 usually have opposite 
parity, independently on the initial condition, both the odd and the even set are 
necessary to completely describe the problem and obtain the correct result. 
In the following paragraphs a compact notation for the space derivatives is intro- 
duced. In order to simplify the reading, the y-derivatives will be indicated with a 
subscript. The temporal derivatives will be indicated explicitly or with a dot. 
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3.2 - Solution to v equation 



3.2.2 Weak formulation and approximate solution to v equation by 
Galerkin method 



Substituting the expansion (3.7) in v equation (3.1) yields 



oo , oo , oo 

= Zl ^c«(^)^".. " ^' Zl ^c„(t)X„ + iaU{y) Y^ c„(t)X„ 

ra=l n=l n=l 



n=l n=l ra=l 

n=l n=l 

The above expression represents an exact form. If only a finite number of modes 
is considered, the equation is not satisfied exactly, so a residual e (dependent on 
the choice of the functions X^) appears at the left hand side 

^ d ^ d ^ 

e{y, t; a, /?) = ^ — c„(t)X„^^ ^ ^^ Z j^c„(t)X„ + iaU{y) ^ c„(t)X„^^ 

n=l ?i=l n=l 



n=l n=l n=l 



2A;2 J^ A;4 J^ 

n=l n=l 



(3.16) 



Galerkin (1915) focused on the problem of minimizing the functional e, so his 
method consists of a variational approach (see also Chandrasekhar, 1961, p. 27- 
32). He showed that the best approximation of the solution is obtained when 
the error is orthogonal to the space of the linearly independent trial functions 
Xn with n = 1,2, ..A^. In this context, given two functions u{y) and v{y) with 
y G f2 = [—1, 1], the following definition of scalar product applies 

{u,v) = u-vdy (3-17) 

Jn 

so, using the above notation, the Galerkin orthogonality condition can be expressed 
as 

(e,X4 = m = l,2,...,iV (3.18) 

Substituting e with its expression (3.16), inverting the integral and the sum signs, 
and taking the time dependent coefficients out of the integral sign, leads to the 

13 



3 - Wave transient analysis: an eigenfunction expansion solution method 

following system of equations 

^ d ^ d ^ 

n=l n=l ra=l 

N N ,2tt( \ 1 ^ 

-iaeJ2cn{t){U{y)X,„Xj-iaJ2cnit){^-^Xr,Xm) - ;^$^c„(t)(X„^^^^,X^) 

n=l n=l n=l 

+ 7^zZ^"W<^"-'^™)-]^zZ^«W<^«'^™) ^' ^ = 1,2,3, ...iV 

n=l ra=l 

(3.19) 

The original partial differential equation is now reduced to a system of A^ ordinary 
differential equations of the first order, where the time dependent coefficients c„(t) 
are the only unknown. The scalar products can be evaluated analytically or com- 
puted by numerical integration and take the following expressions 



2 \2 



A^ A] 



where 



+^YfzftiKln - Knlm) if (u + m) is cveu, n=l^m 
if (n + m) is odd 

-^lll + Kilm iin = m 

^ m,n x^nyyyyj^m) '^n^m,n yo.Z,Z,j 



C0sh{2\n) - C0s{2\n) 

In = . , .^. ^ . .o\ ^ ^™ ^" = ^ 3.23 

Stnh[2An) — Stn[2An) n^co 



For Plane Couette flow, in all the present work the following expression of the base 
flow will be considered 

U{y) = y (3.24) 

so that the other integrals take the expressions 



f/«„=([/(y)X„^^,X„)= (3.25) 



"^^^{Kln - Kalm - 1) " Sy^^^fV^n^m ^^ ("- + "^) ^^ odd 



if (n + m) is even, n ^ m 
if (n + m 
if n = m 

(3.26) 
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3.2 - Solution to v equation 



USn={U{y)Xn,X^)= (3.27) 



if (n + m) is even, n ^ m 

16§^%^ if (n + m) is odd 
if n = -m 



(rUiv') 
USn = (^^^«' Xm) = ^n,m (3.28) 

It is convenient to express the ODEs system (3.19) in a more compact notation: in 
the following, vectors will be indicated either explicitly using braces or with bold 
lower case letters; matrices will be indicated with bold capital letters; constants 
with roman capital letters and physical parameters in italic. The system can be 
written as 

.^ '. ^ A?.n /^.n /-/.n ' 



Re Re Re 



H 



(3.29) 
Hc-Gc = (3.30) 

where D = [-Dm,™] etc., i.e. the element Dm,n is placed at the n'^^ column and at 
the m*^ row of the matrix. H is invertible, so denoting A = H~^ G yields 

c-Ac = (3.31) 

The general solution to the ODEs system (3.31) in the case of matrix M having 
N distinct eigenvalues yUj (either real or complex, see Zill & Cullen, 2005), reads 

c(t) = Ki lie^i* + K2 he^'' + . . . + Kn l^ve^^* (3.32) 

where Ij are the eigenvectors corresponding to /ij and Ki are constants to be deter- 
mined by imposing the initial condition, or alternatively using the Matrix Expo- 
nential notation 

c(t) = e^'{K,} (3.33) 

The coefficients at the initial time, Cq, can be obtained from the inverse transfor- 
mation (3.8) since the initial condition v{t = 0) is known, so finally the solution is 
get by solving the algebraic system 

Co = Ki li + K2 I2 + + Kn In (3.34) 

ho = {K,} = L-' Co (3.35) 

where h(t) = L~^ c(t), and L is the matrix whose columns are the eigenvectors Ij. 
Their linear independence ensures that L is invertible. 
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3.3 Solution to the forced 77 equation 

3.3.1 Choice of a base of orthogonal functions 

Following the same procedure of §3.2.1, we now focus on the normal-vorticity 
equation (3.2), which is forced by the solution -0 of the Orr-Sommerfeld PDE equa- 
tion (3.1). Together with the non-orthogonality of the Orr-Sommerfeld differential 
operator, a resonance phenomenon has been pointed out as one of the reasons for 
large energy transient growths, if there is sufficient wave obliquity (see Gustavsson, 
1991). In order to solve the fj equation a set of normal functions different from 
the one adopted in §3.2.1 is needed, since the second order PDE only requires fj to 
vanish at the boundaries, but not its ffist derivative . The simplest choice for the 
basis functions, here adopted, is the following 

Yn = sin{^ny) n = 1,3,5, ...A^ — 1 {odd set) (3.36) 

Yn = cos{^ny) n = 2,4,6, ...N {even set) (3.37) 

where 

(n + Dtt 
^n=- — ^-^ n= 1,3,5,. ..iV-1 {odd set) (3.38) 

(77/ — 1 JTT 

Cn = ^ — 7^ n = 2,4,6, ...A^ {even set) (3.39) 

Also in this case, note that two sets of eigenfunctions are put together to form 
a unique set, since both are necessary to completely describe the complex- valued 
normal vorticity. The general solution is then obtained as the sum of a particular 
solution r/p and the solution to the corresponding homogeneous equation f]h 

V{y,t)=f]h{y,t)+fjp{y,t) (3.40) 

3.3.2 Weak formulation and approximate solution to fj equation by 
Galerkin method 

Considering the complete equation (3.2), we proceed as done for the normal- 
velocity and expand the solution as follows 

oo 

7){y,t) = J2Ut)Yn{y) (3.41) 

n=l 

Substituting, the equation reads 

oo , oo ^ oo 

J2 ;rMt)Yn + taU{y) J2 Ut)Yn - j^J2 ^" W^«.. 



dt ^^-^ Re 

ra=l n=l n=l 



P ^ .dU{y) 



oo 



(3.42) 



+ ]^ E ^«(^)^" = -'f^^^ E ^n{t)X4y) 

n=l ^ n=l 
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3.3 - Solution to the forced -q equation 



Considering a finite number A^ of terms of the expansion and applying the Galerkin 
method yields 

N . N ^ ^ 

Y, ^^bn{t){Y^, Y„,) + m ^ 6„(t)(f/(i/)F„, Ym) - ^^Y^ K{t){Yn,,, Y„,) 

n=l n=l n=l 

+ -^5^&„(t)(V;,r„) = -2/3 5^c„(t)(^^X„,F„) n,m = 1,2,3, ...,iV 

n=l n=l ^ 

(3.43) 

The scalar products can be evaluated analytically and take the following expres- 
sions 

D*^^^ = {Y^,Y^)=5rn,n (3.44) 

(3.45) 





^m,n ~ \^nyy, J^m) — krfim,n 


Plane Couette flow: 






U*={U{y)Y^,Y^) = 



(3.46) 



if (n + m) is even, or n = m 

n-\-7n-\-l 

^~^U%2fl"^" - if in + m) is odd 

F:n,n = {^^X^,Yr^)= (3.47) 

y— _ m-\-l 

2V2UU-i)^ Z. if ^^ m are odd 

t4 _f4 — if n, m are even 
Introducing the vector and matrix notation, the system (3.43) reads 

1 2P 

h-( - ia U* + — S* - -— D*) b = -i/3 F* c (3.48) 

Re Re ^ — v — ■' 

^ ^ ^ B 

G* 

b-G*b = Bc (3.49) 

which is a non-homogeneuos system of A^ ODEs. The general solution to (3.49) 
consists of the superposition of the solution of the homogeneous system and a 
particular one. Naming $(t) the fundamental matrix of the system ( see Zill & 
Cullen, 2005), then a formal expression for the general solution is 

b(t) = b;,(t) + bp(t) = #(t) k + #(t) f ^-\t) B c(s) ds (3.50) 

Jto 

or, alternatively, using the Matrix Exponential form 

b(t) = Kit) + bp(t) = e^** k + e^'* [ e-^'* B c(s) ds (3.51) 

Jto 
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3 - Wave transient analysis: an eigenfunction expansion solution method 

where k is a vector of constants. The last expressions have the advantage that the 
particular integral vanishes at t = tg, so it is easy to find the constants by imposing 
the initial condition. Unfortunately, a numerical evaluation of the integral can lead 
to non negligible errors, especially for big times, where a product of very large and 
very small terms occurs. In order to make the numerical computation possible, a 
different form for the particular solution is sought. 

Particular solution bp 

Since the solution (3.32) in terms of the expansion coefficients c(t) is a com- 
bination of exponentials and represents the forcing term in (3.49), the following 
particular solution bp is sought 



N 






(3.52) 



where a„j are constants and fij are the eigenvalues of G, through which the forcing 
term is expressed. Yields 



N 



Vp 



{y,t) = J2Kit)Yn{y) 



(3.53) 



ra=l 



Diagonalizing the system (3.31), the coefficients of the normal- velocity result 



' /ime'^i* ' 



c(t) = Lh = L < 



tOif 



ho^e^^' 



ho = L"^ c(t = 0) 



(3.54) 



Substituting the particular solution (3.52) in (3.49) and leads to 






Q. 



fliie^i* + . . . aiNC^"^^ 



ame^^* + ... aNNe''''^ 



( /in, e^i* 1 



BL< 






Ve^^* 



(3.55) 
It is straightforward to find the unknown constants a„j by comparing terms with 
the same exponential factor. This is equivalent to solve the following set of A^ 
algebraic systems 



{^^,I-G*){ 



aij 

«2i 


> = hoj < 




j 


= 1,2, 


..N 


^ajv,_ 




Bh. 









(3.56) 



where I is the identity matrix and B*j are the elements of the matrix BL. As 
usual, the first subscript indicates the row and the second one indicates the column. 
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3.4 - A Matlab® code implementation 



Finally we get the matrix of coefficients column by column as 



aij 

02i 



} = ifijl - G*)-'hoj < 



^2i 



a^j I 






j = l,2,...iV 



(3.57) 



Homogeneous and complete solution b 

The homogeneous solution of (3.49) takes the same form of the of equation 
(3.32). Indicating with /i* and 1* respectively the eigenvalues and eigenvectors of 
the matrix G*, it follows 



b;,(t) = Ci Ite^i* + C2 lle^'i' + . . . + Cn l^e^^* 
Finally the complete solution is 

b(t) = Ci l*e^i* + C2 i;e^2* + . . . + Cn l^.e''^* + \{t) 



(3.58) 



(3.59) 



The unknown constants Cj depends on the initial condition, and can be calculated 
setting t = in the above expression, leading to 



h* = {C,} = L*-i(bo-b 



Op) 



(3.60) 



3.4 A Matlab® code implementation 



3.J1..I Code description 

In order to verify the proposed method and to obtain the numerical solutions, a 
Matlab® code has been developed and used in the further analysis. At present, 
it consists of two scrips: the main code solves the normal-velocity equation and 
calls a function to solve the normal-vorticity equation. Eventually it computes the 
other components of the perturbation velocity, the frequency, the energy growth 
factor and other quantities of interest. 

The main code 

The structure of the main program {main_ivp_ galerkin.m) can be represented 
by the flowchart of Fig. 3.2. Among the other simulation parameters, the number 
of modes can be chosen. The corresponding eigenvalues Aj are computed by a 
separate script through Bisection or Newton-Raphson method, and memorized in 
a .mat ffie which can be loaded by the main program. As seen in §3.2.2, the method 
requires the solution to algebraic systems, so matrix inversions. In particular two 
inversions are needed to obtain the solution v, respectively the one of the matrix 
A and the one of the eigenvectors matrix L. In fact, the ill-conditioning of this 
matrices can influence the accuracy of the computation. In detail, the condition 
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number of matrix A does not represent a problem, while the one of matrix L can 
reach very high values, of the order of 10^^ for some parameters settings as low 
obliquity angles; moreover the condition number increases with increasing A^. This 
fact has not been investigated in details, but is quite similar to that pointed out 
by Schmid & Henningson (2001). The ill-conditioning of L is intrinsic, due to the 
non orthogonality of the Orr-Sommerfold linear operator. In order to guarantee 
a certain level of accuracy, the matrix inversion error is checked and compared to 
a threshold tol. The absolute pointwise error in the solution of a generic system 
A X = b is defined as follows 



err^ 



Ax — b 



(3.61) 




Set parameters: 
k, 0, Re, type, i.e., N, tol 



[Set Base flow] 



Load 
A,: 



TJ(y) dU dHl 



Compute 
eigenfunctions 



7t 5 y ' yy 




Compute 
initial condition 

Co 



Compute 
matrices 




Diagonalize A 



/^i,L 



Compute ODE 
solution 



c(t) 



Compute v{y,t) 




Compute fi{y,t) 



Call 

solve _squire.m 



Postprocessing 



I 




G{t),u{y,t),w{y,t) 
^iy,i),c{y,t) 



Figure 3.2: Structure of the main code main_ivp_galerkin.m. 

The script computes by default the matrix inversion using the Matlab® backslash 
command. Then, the norm (both L2-norm and Loo-norm ) of the error is calculated 
and compared to the tolerance. If the threshold level is exceeded, a different 
method to solve the algebraic system is performed. In particular, the present code 
tries to approximately solve the system using the GMRES method. Both the 
norm of absolute and the relative error are evaluated. Even if a better method 
should be object of future deeper analysis, it has been observed that for every 
parameters combination tried, the minimum relative error is obtained using the 
backslash command, and its order of magnitude spans the range 10^^^ — 10^^°. 

The solution is obtained for every point along the space coordinate y and at all 
time points deflned by the user. Since the time evolution is analytically obtained. 
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3.4 - A Matlab® code implementation 



and since the analytical expressions of the modes is known, the accuracy of the 
method shouldn't depend neither on the space nor on the time discretization. Ac- 
tually, this is not exact: the computation of the coefficients Cq is performed through 
numerical integration (trapezoidal rule), whose order of accuracy is 0{Ay'^). More- 
over, in the case of Plane Poiseuille ffow, the scalar products Um,n, Urn,n and U^ „ 
are evaluated through numerical integration as well, so a sufficient grid spacing is 
required. The computation of the wave frequency or phase velocity can require 
very fine time grids, the motivation will be given in Chapter 4. In general, fine 
grids are necessary whenever a finite difference scheme for derivatives computation 
needs to be applied. In order to ensure an high accuracy of the computation of 
the initial condition, and consequently of the global solution, and at the same time 
allowing the user to obtain the solution only at the desired points, a separate grid 
is used just for the computation of the initial condition's coefficients Cq. 



The fj solving function 

The fj equation is solved by the function solve_ squire. m. The script structure 
is represented by the flowchart of Fig. 3.3, and is quite similar to the one of the 
main program. Here there is no reason to calculate the eigenvalues related to the 
basis functions with another script, being known their exact expression. It can 
be noticed that to obtain the normal vorticity the solution of the A^ algebraic 
systems (3.56) is required, as well as the inversion of the matrix L*. Depending 
on the choice of the simulation parameters, these matrices can be ill-conditioned 
so the same procedure described in the previous section is implemented. During 
the simulation all the error norms are displayed (about the N algebraic systems 
(3.56), only the maximum will be shown). 

Also in the computation of fj, the backslash command has revealed to guarantee 
an high level of accuracy in all simulation performed, being all error norms of the 
order of 10~^^ or lower. 

3.4-2 Rate of Convergence 

In this section, the rate of convergence of the present method to the exact so- 
lutions of the initial value problem will be investigated. The eigenfunction ex- 
pansion method using the Chandrasekhar basis functions was applied to the Orr- 
Sommerfeld modal analysis by Gallagher & Mercer (1962). In that case, the error 
was shown to decrease as 1/A^^ as A^ — > oo, moreover the residue after A^ terms of 
expansion was of order 1/N^. 

Even if the present formulation is different in the fact that the PDEs are reduced 
to systems of ODEs, rather than algebraic equations, it has been verified that the 
convergence ratio of order 1/A^^ as A^ — )■ oo is kept, and it applies at all times, 
for what concerns -0 equation. About fj equation, the method reduces to a Fourier 
series expansion with time-dependent coefficients, and in this case the convergence 
rate is found to be only slightly different from the one of v. 

Both the root mean square (proportional to L2-norm) and the maximum (Loo- 
norm) of the residual are computed. In the following, the residuals as a function 
of the number of modes used for the simulation are reported at two different times 
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Inputs: 
L, /ij, A„ Xn, ho 



Compute 
eigenfunctions 



\* V V 

■^ 1 ■>■ n: ■>■ Uyy 



Compute 
initial condition 



bo 




Compute 
matrices 



Compute [aij] 
constants matrix 




Compute ODE 
particular sol. 



bp(t) 



Diagonalize G 



//* T " 



Compute ODE 
complete sol. 



b(t) 



Compute fi{y^t) 




G*,B,B 




for j = 1 : N 



Figure 3.3: Structure of the function solve _squire.m,. 

and two different values of the obliquity angle, for Plane Couette flow and Plane 
Poiseuille flow. Since the exact solution is not known, the residuals are defined as 
the difference between the solution and an accurate solution computed with 350 
modes. The exact expression is only known at time zero, and the results have been 
confirmed in this case. 



^a{y,t) = \vN{y,t) -VN=350iy,t)\ 
rms{ea){t) 



K, 



N^ 



\ i=l 



max{ea){t) = raeix{e a{y,t)) 



Convergence of v 



(3.62) 

(3.63) 
(3.64) 



The error as function of the number of eigenfunctions is represented in a bilog- 
arithmic plane, so that the slope represent directly the order of the method. It 
can be noticed from Fig. 3.4 and Fig. 3.5 that the error behaves almost like 1/N^ 
as iV — )■ oo for different choices of the parameters and in both norms. Moreover, 
differences between the imaginary and the real part are little. 
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10" 



10" 



- 10" 



10" 



10" 



10" 



Plane Couette flow 




»<- -/m(?)), max[€a) 
— e Re{v), r77is(eo) 



10' 



10^ 



N 



(a) to = 10, Re = 500, k = 6.5, = 30° 



(b) to = 100, Re = 500, fc = 6.5, = 30° 



Figure 3.4: Maximum and mis of the absolute residual of v for Plane Couette flow and symmetri- 
cal (even) initial condition. Continuous line: real part; dashed line: imaginary part. Comparison 
with N-^. 




10" 



Plane Poiseuille flow 
Re = 10^, fc = 2.04 




»<- -/m(?)), t ms {€a)'- 



10' 



10' 



N 



(a) to = 10, Re = 1000, fc = 2.04, = 80° 



(b) to = 100, i?e = 1000, fc = 2.04, 



80° 



Figure 3.5: Maximum and rms of the absolute residual of v for Plane Poiseuille flow and 
symmetrical (even) initial condition. Continuous line: real part; dashed line: imaginary part. 
Comparison with N^^. 
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Convergence of rj 





(a) to = 10, Re = 500, k = 6.5, 



30° 



(b) to = 100, Re = 500, k = 6.5, 



30° 



Figure 3.6: Maximum and rms of the absolute residual of rj for Plane Couette flow and symmetri- 
cal (even) initial condition. Continuous line: real part; dashed line: imaginary part. Comparison 
with N-^. 





(a) to = 10, Re = 1000, k = 2.04, </> = 80° 



(b) to = 100, Re = 1000, k = 2.04, (/. = 80° 



Figure 3.7: Maximum and rms of the absolute residual of ij for Plane Poiseuille flow and 
symmetrical (even) initial condition. Continuous line: real part; dashed line: imaginary part. 
Comparison with N"^. 
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3.4-3 Termwise differentiation and convergence of derivatives 

The analytical expressions of the derivatives of the basis functions are known. 
Thereby, it is easy to investigate the convergence of the series obtained by termwise 
differentiation. In fact, it is known (see Strauss, 1992) that termwise differentia- 
tion is not always possible and if the derivatives of the solution are needed, finite 
differences techniques may be required. The convergence has been estimated ap- 
plying the transform to a known function, i.e. the symmetrical initial condition for 
the normal-velocity. The same profile is used to estimate the convergence for the 
normal- vorticity, even if it won't be used as initial condition in the further analysis; 
thus, in the following figures, the absolute error is indicated as ea^^. 
About the convergence of -0, it can be noticed (see Fig. 3.8) that the series can be 
differentiated termwise three times. The derivative of fourth order converges very 
slowly in L2-norm but the maximum error remains almost constant, which means 
that the convergence is non-uniform (see Fig. 3.10a). The same results have been 
found by Orszag (1971). The derivatives of lower order converge in both norm to 
the exact expression, even if the rate decreases with increasing order of derivation. 
The rate of convergence of v is now slightly less than N^ . 





(a) Error rm.3 



(b) Maximum error 



Figure 3.8: Convergence of v series and its derivatives to the velocity profile defined by ic 



About fj series, it can be noticed from Fig. 3.9 that the termwise differentiation 
can be applied to obtain correct results up to the first derivative. The second 
derivative converges slowly in norm two, but non-uniformly (see Fig. 3.10b). In 
this case the error of fj seems to behave almost like A^~^, so differently from what 
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can be seen in Fig. 3.6 and Fig. 3.7, where the rj solution is forced by the normal- 
velocity V. 




10' 



10 



C3 



10" 



10" 



10 



10 



10 



-e fjyy 



10' 



:o: ' : ' o e — e — 0:0: o: pooo : 




10 



N 




N 

(a) Error rms (b) Maximum error 

Figure 3.9: Convergence of fj series and its derivatives to a profile defined by ic = (1 — y^)^. 




(a) 



30 




(b) 



Figure 3.10: Non-uniform convergence of the four-times termwise differentiated v series and of 
the two-times termwise differentiated fj series to the exact derivatives of zc ~ (1 — y^)^- The 
phenomenon is analogous to the Gibbs one, and subsists at the boundaries y = ±1. 
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Chapter 4 



Wave transient analysis: numerical results 



4.1 Introduction 

In the present section, particular attention is given to the temporal evolution of 
the wave frequency and phase velocity. Recent studies (see Scarsoglio et ai, 2012, 
2009) have been pointing out as through wave frequency investigations precious 
information can be obtained about the different phases that characterize the spatio- 
temporal evolution of a perturbation. 

The importance of a better understanding of the transient live of traveling waves 
relies, among other reasons, in the relation with rapid transition to fluid turbulence. 
In fact, it is believed that the exceptionally large algebraic growth which can occur 
in the disturbance evolution before the asymptotic exponential mode is set, could 
promote a phenomenon known as bypass transition (see e.g. Henningson et ai, 
1994). Actually, it consists of a disturbance growth and breakdown on a timescale 
much shorter than those typical for Tollmien-Schlichting (TS) waves. The term 
is used to emphasize that these scenarios bypass the growth of two-dimensional 
waves and their subsequent secondary instability (see Henningson et ai, 1993). 
It has been shown that the transient algebraic growth can be significant even for 
subcritical values of the Reynolds number, so that finite amplitudes can rapidly 
be achieved, and nonlinear effects can enter into play. We will discuss in the next 
chapter the consequent formation of turbulent spots. 

The frequency temporal evolution has been poorly investigated, probably because 
sheared incompressible flows are viewed as non- dispersive media. Nevertheless, the 
frequency transient behaviour revealed unexpected phenomena, non predictable a 
priori, that being related to the wave phase velocity could have a remarkable 
influence on the main phase speed of a group of waves, in particular on the early 
stage of a natural spot formation, when the non-linear effects can be neglected, 
as shown by Cohen et al. (1991) in the case of boundary-layer flow. Moreover, 
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as pointed out by Kachanov (1994), neither intermittence nor turbulent spots are 
observed in K- or N-regimes of transition when the initial instability wave is strictly 
periodic in time. The "natural" intermittence phenomena and the spot formation 
are usually detected when the perturbations background is more complicated and 
the instability wave has both amplitude and phase modulation in time. 

In this contest, the study of the frequency temporal evolution gains more and 
more importance. From the latest works previously cited emerges that the com- 
plexity of the frequency transient is mainly associated to jumps which appear quite 
far along the temporal history. The normalized time at which the jumps occur have 
been considered as the threshold between the first two phases of a wave life, respec- 
tively the Early transient and the Intermediate transient. The intermediate term 
lasts until the asymptotic exponential energy growth/decay is reached (long term), 
and appears to be the most probable state in a wave life, since on one hand its 
temporal extension is at least one order of magnitude bigger than the early term's 
one and, on the other hand, at the end of this intermediate period the disturbance 
will die or blowup. 

In the following sections, an analysis of the phase speed time evolution will be 
provided; The three phases of a wave life will emerge by taking into account the 
frequency time history of both the normal- velocity and the normal-vorticity; par- 
ticularly, the existence of the intermediate period will be shown. Moreover, a rela- 
tionship between frequency jumps and the achievement of a self-similar asymptotic 
state of the flow velocity and vorticity profiles will be shown. 



4.2 Wave frequency and phase velocity 



4.2.1 Analysis of the v component of flow velocity 

In the present section we take advantage of the solution method proposed in 
§3.2.2 to obtain the time evolution of the frequency and phase speed of the wall- 
normal component of velocity v, varying the three parameters that characterize the 
problem: the Reynolds number, the obliquity angle and the polar wavenumber. 
The high non-stationarity of the phenomenon will emerge, typically a jump is 
observed at a certain time, which is considered to be the threshold between the 
early period and the intermediate one. After this jump the frequency of v is 
characterized by a modulation about a constant mean value for Plane Couette 
flow, for sufficiently high values of the polar wavenumber k. For all the simulations 
performed, A^ = 250 eigenfunctions are used for the solution expansion. 

Numerical computation 

The frequency of the perturbation is defined as the temporal derivative of the 
unwrapped phase 6{y,t;a, P), at a specific spatial point along the y coordinate. 
The wrapped phase 

e^{y, t; a, (3) = arg{v{y, t; a, /?)) (4.1) 
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is a discontinuous function of t in [— vr, +7r], while the unwrapped phase 6 is contin- 
uous and it is obtained by adding multiples of ±27r when absolute jumps greater 
than or equal to tt radians occur. The wave frequency is defined as 



u{y,t;a,f3) 



d0{y,t;a,f3) 



dt 
The phase velocity vector is given then by the dispersion relation 



(4.2) 



c=|k (4.3) 

where k = {cos{(f)),sin{(f))) is the unitary vector defining the polar wavenumber 
direction. The frequency of each signal can be numerically computed at a fixed 
observation point y = ijq. In order to ensure a high accuracy in the results, a fourth 
order centered finite-differences scheme has been used to calculate the first temporal 
derivative. The following scheme applies for the inner points of the defined time 
vector and corresponding phase values {ti,6i): 

ujiU;yo,a,^)=UJ, = ^ = ^^ |^^^ '^' ^ ^ = 3,4, . . . , iV,-2 (4.4) 

where A^^ is the total number of elements of the time and phase vector and At the 
time spacing. It is worth to underline that the proposed method allows the user 
to define arbitrary time and space grids, differently from Runge-Kutta routines, 
where the time step is free to change accordingly to the stiffness of the problem 
and the needed accuracy. This is actually an advantage, because the accuracy 
of the numerical estimated derivatives is affected by the non-uniformity of the 
grid. Setting a uniform time spacing we ensure that the finite-differences scheme is 
actually of the fourth order (see Fertziger & Peric, 1996). Since the scheme stencil 
is made of five points, for the first and the last two points of the vector respectively 
a forward and backward fourth order finite-differences scheme is needed. Indeed, 
the accuracy of the method for these points could be lower. The following schemes 
are applied: 



d^i -25^1 + 48^2 - 36^3 + 16^4 - 3^, 



5 



(4.5) 
(4.6) 

^^*-i = — J^ = TTTT^ (4-7) 






u^i — 


dt 
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Phase velocity temporal evolution 

In the following, the temporal history of the absolute value of the phase velocity 
of the V component, |c(t)|, for the Plane Couette flow is presented, for different 
combinations of the parameters. As one can be notice from figures 4.1, 4.2, 4.3, 4.4, 
the problem, though linear, offers a complex highly non-stationary scenario, which 
is hardly possible to estimate a priori. Two different periods for the phase velocity 
(and so the frequency) temporal evolution can be observed, the Early term and 
the Intermediate term. The last one ends when the perturbation energy growth 
factor reaches the exponential asymptotic trend, and will be investigated in detail 
in the following section. The transition between the early and the intermediate 
transient appears to happen in a narrow time window, and it is often characterized 
by an abrupt jump to a higher mean value which is maintained throughout the 
rest of the perturbation's life, indicated as c or u for the phase velocity and the 
wave frequency respectively. 

Another important observation is that the asymptotic value generally is not 
a constant one for PCf, but a modulation characterized by a specific period T^ is 
present. For this reason it is suitable to refer to frequency asymptotic mean values. 
A motivation of this fact will be provided in the next paragraph. As it clearly 
appears from Fig. 4.1b and Fig. 4.2b, the oscillations amplitude can be significant 
with respect to the mean value, specially for low wavenumbers. With increasing 
k, the frequency temporal evolution appears to shift from peaks-characterized to 
sinusoidal. 

Moreover, it has been verified that exists a certain threshold in the wavenumber, 
denoted with kj, below which neither jump occurs nor phase velocity modulation is 
observed, but the wave frequency experiences a monotonic decay to the zero value, 
after a transient evolution (see Fig. 4.1a, 4.2a and Tab. 4.3). 

In Tab. 4.1 and Tab. 4.2 values of nondimensional time at which the jump occurs are 
reported. Since the jump is spread within a certain time window, a strict definition 
of Tj is not provided, so the normalized time corresponding to the frequency peak 
that typically occurs after the jump is considered, in the present work, as an 
index of the end of the early transient. Even if the temporal evolution of the 
wave frequency varies with respect to the simulation parameters and the "jump" 
itself shows different shapes, some general trends can be observed. In fact, Tj 
seems to increase with increasing Reynolds number (Fig. 4.3a) and obliquity angle 
(Fig. 4.4a), while it decreases with increasing polar wavenumber (Fig. 4.3b). 
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Figure 4.1: Temporal evolution of the absolute value of the phase velocity, calculated from v, 
for PCf at Re = 500, = 0° and sym. initial condition, defined by Vq = (1 — y^)^. The polar 
wavenumber covers the range k G [0.05,25], uniformly distributed in the logarithmic space. The 
fixed observation point is j/o = 0-5. 
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Figure 4.2: Temporal evolution of the absolute value of the phase velocity, calculated from v, 
for PCf at Re = 500, (p = 45° and sym. initial condition. The polar wavenumber covers the 
range k G [0.05,25], uniformly distributed in the logarithmic space. The fixed observation point 
is yn = 0.5. 
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Figure 4.3: Temporal evolution of the absolute value of the phase velocity, calculated from v, 
for PCf for k — 10, (j) — 45° and sym. initial condition. The Reynolds number covers the range 
Re e [10,10000], uniformly distributed in the logarithmic space. The fixed observation point is 
j/o = 0.5. (a) Detail of the Early transient and frequency jumps; (b) Intermediate transient. 
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Figure 4.4: Temporal evolution of the absolute value of the phase velocity, calculated from v, 
for PCf for k = 6.5, Re = 500 and sym. initial condition. The obliquity perturbation angle 
covers the range ip e [0°,90°], uniformly distributed. The fixed observation point is j/o — 0.5. 
The case (f> — 90° is not represented since the orthogonal wave is stationary for all k. Re and 
initial condition. 
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Table 4.1: Frequency jump nondimensional time Tj for various combination of the simulation 
parameters, for Re = 500 and sym. initial condition. Tj is considered as the time at which the 
frequency maximum value, typically located just after the jump, occurs. 
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Table 4.2: Frequency jump nondimensional time Tj for various combination of the simulation 
parameters, for Re — 5000 and sym. initial condition. 



Intermediate Term and Long Term behaviour 

An interesting aspect of the nonmodal analysis, hardly ever considered in the 
past, is the possibihty to investigate on how the asymptotic state is reached. In 
the present and in the following sections some new results will be presented. The 
trends of the mean values of the asymptotic phase velocity for various combinations 
of the simulation parameters are reported in Fig. 4.5a and Fig. 4.5b. Accordingly 
to known results from the cited literature, the phase velocity asymptotic mean 
value increases with increasing k, Re and decreasing 0. The asymptotic value 
doesn't depend neither on the initial condition nor on the observation point yo, 
but only on the spectrum of v for that particular parameters combination. In fact, 
the frequency corresponds to the real part of the least damped eigenvalue of the 
Orr-Sommerfeld operator, while the damping is given by the imaginary part. Even 
if the present method is developed to study the temporal evolution of perturbations 
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from an initial value problem, the spectra af v can easily be obtained by computing 
the eigenvalues of the matrix A (see §3.3.2). 

About the mean values of the asymptotic frequency u, trends are shown in Fig. 4.6a 
and Fig. 4.6b. It can be noticed that for high values of k, the general tendency can 
be approximated by the relation u = kcos{(f)). This trend was observed for Plane 
Poiseuile flow and for Wake flow, as well, by Scarsoglio et al. (2012). 
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Figure 4.5: (a) Trends of the absolute mean value of phase velocity, calculated from v, for 
PCf for Re = 500, (j> = {0°,45°,80°} and sym. initial condition. The polar wavenumbers are 
uniformly distributed in the logarithmic space. The fixed observation point is yo = 0.5. (b) 
Absolute mean values of c, calculated for PCf for (p = 45°, Re — {50, 500, 5000} and sym. initial 
condition. These asymptotic results are independent on both the initial condition and yg- 



34 



4.2 - Wave frequency and phase velocity 



120r 



100 



13 60 



0° Plane Couette flow 

45° fle = 500, yo=0.5, Symi.c. 

80° Intermediate and far period 

\^ ^ tjump I 




100r 



90 



80- 



-Re 
■Re 
■Re 



50 

500 

5000 



Plane Couette flow 
= 45°, 1/0 = 0.5, Sym 



Intermediate 



md f ar period 



120 




(a) 



(b) 



Figure 4.6: (a) Trends of the absolute mean value of the frequency, calculated from f), for PCf for 
Re — 500; (f) — {0°,45°,80°} and sym. initial condition. The polar wavenumbers are uniformly 
distributed in the logarithmic space. The fixed observation point is j/q = 0.5. (b) Absolute mean 
values of Co, calculated for PCf for — 45°, Re ~ {50,500,5000} and sym. initial condition. 
These asymptotic results are independent on both the initial condition and y^. 



As seen in the previous section, there exists a threshold for k below which the 
temporal evolution is characterized by a frequency decay to zero, after an early 
transient. To better observe these conditions, in Fig. 4.7 the phase velocity evo- 
lution has been traced on a semilogarithmic plane. Also for these cases there is 
no influence of the initial condition on the far periods. Moreover it can be noticed 
that for k < kj there is a phase of exponential decay leading to a stationary state. 
The phenomenon has been observed by Gallagher & Mercer (1962) through a modal 
analysis (2D). They discovered that below a certain value of aRe all the eigenval- 
ues were real. Increasing Re, a threshold level is reached, where the least damped 
eigenvalues are real and coincident; they split into a complex conjugate pair for 
larger values of the Reynolds number. The authors pointed out the abruptness of 
the transition, as well. Here various simulations have been performed in order to 
verify the precision of the method and of the numerical code for the asymptotic 
solution computation. The results have been compared to those of the cited au- 
thors, and an excellent agreement has been found. In addition, a generalization 
for the three-dimensional case is shown (see Fig. 4.8b). 

Another comparison have been done with the results obtained by Orszag (1971), 
who developed a method for the modal analysis based on a Chebyshev polynomials 
expansion. For Plane Poiseuille flow at Re = 10000, a = 1 and = 0° he found the 
following value for the unstable eigenvalue: c = 0.23752649 + 0.00373967i. With 
N=250 we obtain the value of c = 0.23752629 + 0.00373964? (with Ay = 10-^ 



120 
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N = 250). The same accuracy is assured for the the other eigenvalues. 
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Figure 4.7: Temporal evolution of the absolute value of the phase velocity, calculated from v, 
for PCf for (j) = 45; Re = 1000. An abrupt variation in the temporal evolution trend can be 
observed once k — 0.105 is exceeded. Dot-dashed line: asym initial condition; continuous line: 
sym i.e. 
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Table 4.3: Values of threshold wavenumber kj below which the asymptotic frequency tends to 
zero, i.e the wave tends to a stationary (damped) state at high times. Here kj is reported for 
three different values of obliquity angle and for ten values of Reynolds number, spanning four 
decades, uniformly distributed in the logarithmic space. 
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Figure 4.8: Trends of the threshold wavenumber kj for PCf, as a function of Re for = 
{0°,45°,80°} (a). Comparison with results found in literature (Gallagher & Mercer, 1962) for 
the bidiniensional case (0 = 0). It is worth to notice that here the same trend is found for the 
three-dimensional cases. 
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Figure 4.9: Trends of the phase velocity asymptotic period T^ for PCf, as a function of k for 
(f> = {0°,45°,80°} and Re = 500 (a). Re = {50,500,5000} and = 45° (b). The observation point 
yo and the initial condition does not influence the period of the modulation. 
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We find interesting to investigate the trend and the reasons of the phase velocity 
asymptotic modulation, since no detailed literature is found about this topic. The 
modulation is characterized by a period T^ that decreases with increasing polar 
wavenumber, according to an exponential law, for sufficiently high values of k 
(Fig. 4.9a). Tc increases with increasing obliquity angle as well, while the influence 
of Re is weak at high k. 

Even if further investigations are needed to verify these results, we underline the 
total agreement with the results obtained by direct numerical integration of (3.1) 
and (3.2) with the methods of lines and Runge-Kutta ODE solver. In the following 
we try to provide a motivation about their non-contradictory nature with the modal 
theory. The spectrum of channel flows (for sufficiently high values of k) is generally 
composed by three branches, whose label A,P,S were given by Mack (1976). In 
the case of Plane Poiseuille flow all the three branches are present and correspond 
respectively to wall modes (c^ — )■ 0), center modes (c,. — )■ 1) and highly damped 
modes (cr — )■ 2/3). In the case of Plane Couette flow, the spectrum does not 
contain a P branch but it has two A branches, composed by complex conjugate 
eigenvalues. 

The solution of the initial value problem generally contains the contribute of all 
the frequency components, as can be clearly seen from the solution of the ODE- 
reduced velocity equation (3.32). Considering the spectrum of PCf, we observe 
that there are two least damped, complex conjugate, eigenvalues (with the same 
damping rate and opposite real frequency). It has been observed, as well, that 
the coefficients ho are complex conjugate; the same does not apply, however, to 
the final solution where the coefficients are mixed up by the matrix L and even- 
tually multiplied by the corresponding shape function X„. In this case the phase 
of the final solution results oscillating and so the frequency. The mean value of 
the asymptotic frequency corresponds exactly to the real part of the least damped 
eigenvalues pair. 

To support this motivation, simulations have been performed for PPf, for param- 
eters combinations where the least damped eigenvalue is unique. As expected, 
for these cases no phase velocity oscillations are observed. In the following, some 
spectra examples for both Plane Couette flow and Plane Poiseuille flow are pro- 
vided. We remind that in the classical modal analysis the solution, e.g. the normal 
velocity, is expressed as 

v{y,t)=d{y)e-"'' (4.9) 

where u = Ur + ioOi is the complex eigenvalue, so the real part Ur represent the 
wave frequency, while the imaginary component Ui is damping rate. 
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According to the convention adopted in the present work (see the exponential 
terms in (3.32)), the damping rates of the single modes are given by the real part 
of the eigenvalues, /x^ = oji and the frequencies correspond to the imaginary part 
with changed sign /Xj = —u^- In order not to be confusing, we will express the 
spectra with the classical convention found in literature, in terms of phase velocity. 
Only the least damped eigenvalues are represented in Fig. 4.10 and Fig. 4.11. 
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Figure 4.10: (a) Spectrum of both the Orr-Sommerfeld and the Squire operators for PCf at 
Re — 1000, (p = 45° and k = y/2. The least damped eigenvalue, respectively for v and for i), is 
explicitly reported, (b) Spectra for PCf at Re = 1000, = 45° and k = 0.08. Here the polar 
wavenumber is below the critical value, in fact we observe that the least damped eigenvalue is 
real. This means that v tends to a stationary damped state as t -^ oo. 

As one can notice, the spectra of the Orr-Somerfeld (-0) and Squire (f/) operators 
are usually different, and for Plane Couette flow in many cases the least damped 
eigenvalue belong to the set of the Squire operator (see for example Fig. 4.10a). 
this fact is found to have an influence on the dynamic of the system which has not 
been taken in account yet, as shown in the next section. 
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Figure 4.11: (a) Spectrum of both the Orr-Sommerfeld and the Squire operators for PPf at 
Re ~ 1000, 4> = 45° and fc = 1. This is an instable configuration: an eigenvalue with positive 
imaginary part exists, (b) Spectra for PCf at Re = 1000, (j) — 45° and k = 10. With increasing 
fc, or Re, the eigenvalues exact computation becomes difficult. The issue concerns the sensitivity 
of the spectrum to small perturbations (e.g. the computer finite precision). This is a property 
of the linear operator rather than a property of the numerical scheme, and the junction point of 
the three branches exhibit the largest sensitivity (see Schmid & Henningson, 2001). However, the 
validity of the solution is not compromised, since the problem doesn't affect the first eigenvalues. 
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4-2.2 Behaviour of the vorticity component fj and global considera- 
tions 

In order to understand the behavior of the complete solution, the normal vor- 
ticity must be considered or, alternatively, the other components of perturbation 
velocity u and w. In the present section the phase velocity of the vorticity signal 
is investigated and some new results are presented. The same fourth order finite- 
difference scheme introduced in §4.2.1 is used for the computation of the phase 
first derivative. 

Temporal evolution of the fj phase velocity 

As seen in the previous section, the non-modal analysis allows to observe the 
complete life of a perturbation, from the early transient to the asymptotic state 
predicted by the modal analysis. However, in the past the same frequency for both 
the normal velocity and the normal vorticity (or, similarly, for the three components 
of velocity) was usually considered by the authors dedicated to the modal analysis. 
As pointed out by Schmid & Henningson (2001), only the particular solution of the 
modal Squire equation has the same frequency of v. In the following analysis, the 
role of the homogeneous part fjh in the frequency temporal evolution of fj is shown. 
Moreover in the following section we will focus on the evolution of the velocity and 
vorticity profiles along the y coordinate, and their correlation with the frequency 
time history. 

As shown in Fig. 4.12a, for sufficiently high values of Re, the phase velocity of fj 
(c^, in the following) has approximately the same evolution of the phase velocity 
of V analized in the previous section, c„ in the following, even if it is clear that they 
are not coincident. To be more precise, a certain lag in the first jump time Tj is 
observed. 

Figure 4.12b reveals an interesting aspect: the frequency of fj experiences a second 
jump, after which it reaches the asymptote predicted by the modal theory. In the 
time window between the these two jumps the mean value of c^ is about the one of 
Cv, and this phase of the wave life can also last several time units, depending on the 
parameters. In fact, the time at which the second jump, Tj2, occurs increases with 
increasing Re and with decreasing k, while about the infiuence of the obliquity 
angle, we observe that Tj2 increases with increasing for high k, but the opposite 
trend occurs at lower wavenumbers (see Tab. 4.4). 
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Figure 4.12: Temporal evolution of the absolute value of the phase velocity, calculated from fj 
(continuous line), for PCf for /c = 10, = 45°, Re — {27,193,1389} and sym. initial condition. The 
fixed observation point is y^ — 0.5. (a) Detail of the early transient and first frequency jump; 
(b) Intermediate and far transient: the second jump can be observed. The phase velocity of v is 
shown with dot-dashed lines. 
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Table 4.4: Frequency jump nondimensional time Tj2 for various combination of the simulation 
parameters, for Re = 500 and sym. initial condition. Since the transition to the asymptotic value 
of Crj can be more or less smooth, Tj2 is considered as the time at which the frequency peak, 
typically located just after the jump, occurs. 
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Intermediate Term and Long Term behaviour 

In order to understand the physical reasons why Cr, experiences two jumps during 
its temporal evolution, we take advantage of the mathematical formulation intro- 
duced in §3.2 and §3.3. The Squire equation (3.2) is forced by the solution -0 of the 
Orr-Sommerfeld PDE (3.1). The general solution can be expressed as fj = 7]^ + fjp 
as shown in §3.3.2. The particular solution contains the same eigenvalues spectrum 
of the forcing v, while the spectrum of the Squire operator (the homogeneous part) 
fih is different. Thinking about a generic forced linear system, it is clear that the 
asymptotic solution has the same frequency of the forcing term if its amplitude is 
constant. If the forcing term itself is damped, the asymptotic frequency depends 
on the damping of both the forcing term and the homogeneous solution. If the 
damping rate of the forcing term is higher than the one of the homogeneous oper- 
ator, the frequency for t — t- oo will be the "natural pulsation" of the system. 
Here the system is far more complicated but the same phenomenon is observed; 
for several configurations of the parameters, looking at the spectra (e.g. Fig. 4.10, 
Fig. 4.11) one can notice that the least damped eigenvalue belongs to the Squire 
set. In these cases, the second jump of c^ occurs. Tj2 depends on the initial coef- 
ficients of the series ho (i.e. on the initial condition) and on the ratio of the real 
part of the eigenvalues /x, to /i*, and can be qualitatively considered as the end of 
the intermediate term and the beginning of the asymptote, as can be seen from 
the trends of the kinetic energy growth rate in Fig. 4.13. 
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Figure 4.13: Temporal evolution of the kinetic energy growth rate (blue line ) and its derivative 
(green line) for PCf with k = 6.5, (p = 45°, Re = 500 and sym. initial condition. It is evident a 
correlation with the frequency jumps shown in the following plots. 
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Figure 4.14: (a) Temporal evolution of the absolute value of the phase velocity of 77, fjh, fjp and 
comparison with c„ for PCf with k — 6.5, (f> = 45°, Re = 500 and sym. initial condition (upper 
plot). In the lower plot the spectrum for the same configuration is shown, (b) Phase velocity for 
Plane Poiseuille flow, same parameters configuration. In both cases the least damped eigenvalue 
belongs to the Squire set. 
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For Plane Couette flow, the same modulation observed in the phase velocity of 
V is found in the fj component, even if the characteristic amplitude and the period 
are generally different. The same motivation discussed in the previous section 
applies, the eigenvalues of fj for this type of flow are complex conjugate, indeed. 
The frequency of this modulation appears to be generally higher than the one of 
Ct, , supporting the fact that this modulation is related to the imaginary part of 
the least damped /ij (for c^) or /i* (for c,,), as can be inferred by looking at the 
spectra, since usually |3fJ(c^,)| < |3fJ(c^)| for the least damped. The trend of the 
asymptotic frequency is reported in Fig. 4.15a and Fig. 4.15b, where one can see 
that as A; — 7- oo the difference between c^ and c^ tends to vanish. Anyway the 
general trend of the two frequencies, varying the parameters, is approximately the 
same. 
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Figure 4.15: (a) Asymptotic absolute values of c^, for PCf with Re = 500, (/> = {10°, 45°, 80°}. 
The polar wavenumbers are uniformly distributed in the logarithmic space, (b) Asymptotic 
absolute values of w^. Comparison with the trends of c^,, Cjy, plotted with dot-dashed line. 
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4.3 Velocity and vorticity profiles, similarity considerations 
and solutions in the physical space 

4.3.1 Profiles of V, 17 and their similarity properties 

In this section, the temporal evolution of the normal vorticity and velocity pro- 
files along the y coordinate is investigated. We observe that the frequency jumps 
previously introduced are strictly related to the spatial distribution of the solu- 
tions V and fj, i.e. to the distribution of the complete flow fleld, in the wavenumber 
space. In flgures 4.16-4.22 the velocity and vorticity proflles are reported, together 
with the phase velocity time history and the evolution of the flrst derivative of the 
kinetic energy growth rate. Actually, the quantities analyzed in the following are 
the modules of the complex-valued solutions v and fj 

\v\ = yM^dT^H) \fi\ = y^Wfj + '^^ (4.10) 

The module of the general quantity in the wavenumber space can be related to the 
solution in the physical space. In fact, taking advantage of linearity, the inverse 
transform for a single wave reads (Criminale, 2003) 

v{x,y,z,t) = ^ [v{y,t)e'^^+'^' + v*{y,t)e-'''^-'^^] (4.11) 

r/(x, y, z,t) = ^ [v{y, t)e^"^+^^" + fj* {y, t)e-^"^-^^^] (4. 12) 

where the * sign represents the complex conjugate. Hence, the sum of the flrst 
complex quantity at right hand side and its conjugate represents the real distur- 
bance quantity in the physical space; the same applies for u and w, derived from 
(2.25) and (2.26). Since the complex conjugate values can be easily obtained once 
V and fj are computed, this is a convenient way to express the solution. 
The explicit relation between the real and imaginary part of the solutions and the 
quantities in the physical space is derived from the above expressions 

v{x, y, z, t) = ^i) cos{ax + f3z) — ^i) sin{ax + f3z) (4-13) 

fi{x, y, z, t) = ^fj cos{ax + f3z) — ^fj sin{ax + f3z) (4-14) 

The proflle along the coordinate y of the module \v\{y,t) or \fi\{y,t) indicates the 
envelope of the maxima of v, or fj, at a flxed point {x, z) 

\v\{y,to) = raax{v{x,y,z,to)} (4.15) 

x,z 

\v\iy,to) = max{fi{x,y,z,tQ)} (4.16) 

x,z 

The following flgures show how the temporal evolution of the disturbance phase 
velocity is closely related to the spatial distribution. To be more precise, the 
solution in terms of modules seems to achieve a self-similarity in time, when the 
frequency becomes constant. In fact, in these conditions the proflles coincide if 
normalized with their Loo-norm (the maximum along y) or, similarly, with the 
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L2-norm. This means that the space-dependent and time-dependent parts of the 
solution are separable. 

\v\(v f) 

' '^^' ^ - fiy) Self - similarity (4.17) 



(raa,Xyi)){t) 

Usually the component of normal velocity is found to achieve this condition after 
Tj, the time at which the first frequency jump occurs, as shown in Fig. 4.16-4.17 
for Plane Couette flow, and Fig. 4.18-4.21 for Plane Poiseuille flow. The vorticity 
profile continues to evolve until the second phase velocity transition occurs, for 
t = Tj2. For the cited cases, we observe that for PCf the {fjl spatial distribution 
varies quite smoothly (Fig. 4.17) while for PPf an abrupt variation in the parity of 
the profile occurs (the double hump of the modules correspond to odd profiles in 
the physical plane), as shown in Fig. 4.19. 

An interesting case is shown in Fig. 4.20-4.21, where the sudden profile change, 
and the associated second frequency jump are experienced by the velocity compo- 
nent rather than the vorticity one. This is probably due to the infiuence of the an- 
tisymmetrical initial condition on the early and intermediate wave transient. This 
influence may be related to the parity of the asymptotic state, which is indepen- 
dent on the initial condition. It is also interesting to notice that the intermediate 
phase, starting after the flrst jump, is usually very close to similarity conditions; in 
this term, the phase velocities of the two signals are nearly coincident. Moreover, 
we underline that the intermediate transient is, in addition to the early period, 
the most relevant term in a perturbation's life. Indeed, in the introduced cases Tj2 
occurs when the wave kinetic energy is extremely small (the last jump represents 
the beginning of the asymptotic conditions). 

A connection between the periodic frequency modulation observed for Plane 
Couette flow in §4.2.1 and the spatial distribution is pointed out in Fig. 4.22; 
here it should be noticed that a periodic continuous variation in the (normalized) 
proflles of |-0| and \fj\ occurs. Actually, the periodic change happens in the channel 
central region, while the near-wall region remains unchanged and self-similar. The 
corresponding case in the physical space is reported in Fig. 4.30. 
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FiGURE 4.16: Plane Couette flow early transient for Re — 500, (p = 45°, k ^ 6.5 and asym. initial 
condition. Top: phase velocity temporal evolution for the v and f] disturbance (respectively, blue 
continuous line and blue dashed line) and first derivative of the kinetic energy growth rate (green 
line). Middle: profiles of the modulus of v, normalized with respect to the maximum (Loo-norm). 
Bottom: profiles of |?7|/||?7||oo- 
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Figure 4.17: Plane Couette flow intermediate and far term for Re — 500, — 45°, k — 6.5 
and asym. initial condition. Top: phase velocity temporal evolution for the v and i) disturbance 
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growth rate (green line). Middle: profiles of the modulus of v, normalized with respect to the 
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Figure 4.18: Plane Poiseuille flow early transient for Re ~ 500, (p — 45°, k — 6.5 and sym. initial 
condition. Top: phase velocity temporal evolution for the v and fj disturbance (respectively, blue 
continuous line and blue dashed line) and first derivative of the kinetic energy growth rate (green 
line). Middle: profiles of the modulus of v, normalized with respect to the maximum (Loo-norm). 
Bottom: profiles of |?7|/||?7||oo- 
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Figure 4.19: Plane Poiseuille flow intermediate and far term for Re — 500, </> = 45°, k = 6.5 
and sym. initial condition. Top: phase velocity temporal evolution for the v and i) disturbance 
(respectively, blue continuous line and blue dashed line) and first derivative of the kinetic energy 
growth rate (green line). Middle: profiles of the modulus of v, normalized with respect to the 
maximum (Loo-norm). Bottom: profiles of |?7|/||?7|loo- It is interesting to notice the abrupt 
transition to the final state of the vorticity profile, which occurs in correspondence to the second 
frequency jump. Both the states before and after Tj2 seem to have similarity properties. 
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Figure 4.20: Plane Poiseuille flow early transient, for Re = 500, cj) — 45°, k ~ 6.5 and asym. 
initial condition. Top: phase velocity temporal evolution for the ij and r) disturbance (respectively, 
blue continuous line and blue dashed line) and first derivative of the kinetic energy growth rate 
(green line). Middle: profiles of the modulus of v, normalized with respect to the maximum 
(ioo-norm) . Bottom: profiles of |?7|/||??||oo- 



52 



4.3 - Velocity and vorticity profiles, similarity considerations and solutions in the physical space 



o 

s 
o 



0.75 



0.72 



0.69 



0.66 



0.63 



A 


\ 




1 


1 1 


1 1 

Plane P oiseuille flow 
k = 6.5, = 45° 
Re — 50Q, Asytn i.e. 


- 


\ 

/ 
/ 


\ 

/ 


^^J^^ 











_^ 


Cjj 




/ V 


^_^ 


1 
1 
1 1 










Intermediate / Long term 



x10" 



2.5 



-to 



-2.5 



50 



100 



150 



200 



250 



300 



5i 








100 



150 



200 



250 



300 



0.5 



?3i - 



-0.5 








50 



100 



150 
f. 



200 



250 



300 



Figure 4.21: Plane Poiseuille flow intermediate and far term for Re — 500, (f) = 45°, k — 6.5 
and asym. initial condition. Top: phase velocity temporal evolution for the v and fj disturbance 
(respectively, blue continuous line and blue dashed line) and flrst derivative of the kinetic energy 
growth rate (green line). Middle: profiles of the modulus of v, normalized with respect to the 
maximum (Loo-norm). Bottom: profiles of |?)|/||?)||oo- The difference from the case of Fig. 4.19 
should be noticed: indeed, here \v\ stabilizes to the final self-similar state after the vorticity 
component. This state is announced by a the second transition of c„ rather than the "usual" 
one of Cjj. This is likely due to the combination of the symmetry properties of both the initial 
condition (here antisymmetrical) and the asymptotic states. 
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Figure 4.22: Plane Couette flow intermediate and far term for Re = 500, (p = 45°, k = 0.63 
and sym. initial condition. Top: phase velocity temporal evolution for the v and fj disturbance 
(respectively, blue continuous line and blue dashed line) and first derivative of the kinetic energy 
growth rate (green line). Middle: profiles of the modulus of v, normalized with respect to the 
maximum (Loo-norm). Bottom: profiles of |?y|/||?y|loo- This is an interesting case characterized 
by large frequency oscillations; the modulation of phase velocity is related to a periodic variation 
of the normalized profiles of both \v\ and \fj\. 
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4.3.2 Maxima of kinetic energy for Plane Couette flow 

Even if the focus of the chapter is on the wave frequency and the similarity 
properties of the velocity profiles, it is thought to be appropriate to include this 
little paragraph about the maxima gained by the kinetic energy during the per- 
turbation's life. The maps of figures 4.23-4.26, together with the evolution of 
the real normalized velocity and vorticity fields introduced in the next paragraph, 
contribute to gain understanding of the complete scenario. 

It is known that in the early and intermediate terms even large transient growths 
can be experienced by the components of flow velocity, vorticity, and by the kinetic 
energy. The normalized kinetic energy density G defined in §2.1.3 can effectively 
measure the transient growth for a perturbation with prescribed initial condition. 
Following the definition by Criminale (2003), an asymptotically stable configura- 
tion is called algebraically unstable if G > for some t > 0; algebraically stable if 
G < for all time; algebraically neutral if G = for all time. The reasons for the 
algebraic growth are mainly three. First, the non-orthogonality of the eigenfunc- 
tions, as shown by Schmid & Henningson (2001). Secondly, a possible resonance 
between the Orr-Sommerfeld and the Squire damped exponential modes can occur, 
as shown by Benney & Gustavsson (1981). However, the resonance does not occur 
for the boundary layer. The last reason deals with the presence of a continuous 
spectrum (so, it only applies to unbounded flows), see the work by Criminale & 
Drazin (1990). 

In the following, the maxima of G are traced as a function of the obliquity angle. 
The nondimensional time at which the maxima occurs is reported as well. Curves 
for six values of Reynolds number (Fig. 4.23 and Fig. 4.25) and polar wavenum- 
ber (Fig. 4.24 and Fig. 4.26) are shown, for both the symmetrical (Fig. 4.23 and 
Fig. 4.24) and the antisymmetrical initial condition (Fig. 4.25 and Fig. 4.26). It 
is interesting to notice that, for fixed Re and /c , it is not generally true that the 
maximum occurs for = 90°. This is evident from Fig. 4.26. 
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Figure 4.23: Upper plot: maxima of the kinetic energy density G as a function of the obhquity 
angle (/>, parametrized with Re. for Plane Couette flow with fc = 6.5 and antisymmetrical initial 
condition. Lower plot: nondimensional times corresponding the the maxima of G. The values of 
Reynolds number are uniformly distributed in the logarithmic space. 
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Figure 4.24: Upper plot: maxima of the kinetic energy density G as a function of the obhquity 
angle 4>, parmetrized with fc, for Plane Couette flow with Re — 500 and symmetrical initial 
condition. Lower plot: nondimensional times corresponding the the maxima of G. The values of 
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Figure 4.25: Upper plot: maxima of the kinetic energy density G as a function of the obhquity 
angle (^, parametrized with Re. for Plane Couette flow with fc = 6.5 and antisymmetrical initial 
condition. Lower plot: nondimensional times corresponding the the maxima of G. The values of 
Reynolds number are uniformly distributed in the logarithmic space. 
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Figure 4.26: Upper plot: maxima of the kinetic energy density G as a function of the obhquity 
angle (^, parmetrized with /c, for Plane Couette flow with Re = 500 and antisymmetrical initial 
condition. Lower plot: nondimensional times corresponding the the maxima of G. The values of 
polar wavenumber are uniformly distributed in the logarithmic space. 
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4.3.3 Wave solutions in the physical space 

In the following, some solutions among those introduced in §4.3.1 are inverse- 
transformed using the relations (4.11) and (4.12) to obtain the quantities in the 
physical space. A few visualizations in the xy plane are here reported (figures 4.27 
to 4.30), for the same cases introduced at the end of the previous section, in order 
to clarify the physical meaning of the module of the complex quantities in the 
wavenumber space, and to observe the behavior of the flow quantities in the real 
three-dimensional space. For all the following flow visualizations a variable color 
scale is adopted to represent at all times the solution, that consequently has to be 
intended as normalized to its maximum value. 
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Figure 4.27: Visualization of the wall-normal perturbation velocity v{x,y, ZQ,tQ), xy plane, for 
Plane Couette flow with the same parameters of Fig. 4.16. It should be noticed that the x- 
component of the phase velocity takes the same sign of the base flow. Remind that for the 
adopted conventions the base flow is oriented as the longitudinal axis x; with reference to the 
figure, the upper wall moves to the right, the lower to the left. Red: maximum (positive); blue: 
minimum (negative). 
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Figure 4.28: Visualization of the wall-normal perturbation vorticity 77(2;, y, zq, to), xy plane, for 
Plane Poiseuille flow. The same parameters of Fig. 4.20 are set, to allow a comparison. From these 
slices, the disturbance seems to move in the streamwise direction, meaning that the a;-component 
of phase velocity has the same sign of the base flow, as observed in Couette flow. Remind that 
the disturbance direction is deflned by the polar wavenumber vector. Red: maximum (positive); 
blue: minimum (negative). 
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Figure 4.29: Visualization of the wall-normal perturbation vorticity fi{x,y, ZQ^to), xy plane, for 
Plane Poiseuille flow. This case correspond to Fig. 4.21, and the effects of the transition to the 
asymptotic conditions, through the second frequency jump of w,,, is shown. Red: maximum 
(positive); blue: minimum (negative). 
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Figure 4.30: Visualization of the wall-normal perturbation velocity v{x,y, zo,to), xy plane, for 
Plane Couette flow. This low-wavenuniber condition was introduced and analyzed in Fig. 4.22. 
From the physical space, the high non-stationarity of the whole scenario is evident. The two 
opposite moving layers, where the largest disturbances reside, appears to experience a periodic 
attachment and consecutive detachment. The images sequence represents one period in the 
modulation of Cy, the maximum of phase velocity occurs in the conditions (a) and (d). Red: 
maximum (positive); blue: minimum (negative). 
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Chapter 5 



Wave packets linear evolution 



5.1 Introduction 



In the present chapter the evolution of hnear wave packets is investigated. The 
aim of this study, as stated in the Introduction, is to emphasize the role of the 
linear mechanisms in a scenario preceding the breakdown and the transition to 
turbulence. To be more precise, the focus will be on the Plane Couette flow, 
extensively studied in the previous chapters, and on Blasius boundary layer flow 
(Bbl, in the following). The wave solutions for the latter are obtained with a 
Runge-Kutta code by numerical integration of the Orr-Sommerfeld and Squire 
PDE equations by the method of lines (see Ames, 1977). 

Bypass transition and turbulent spots 

Although the Plane Couette flow is stable to infinitesimal perturbations for all 
values of the Reynolds number, experimental evidences showed that for sufficiently 
high values of Re the fiow becomes turbulent. This process is observed in bounded 
fiows and in Bbl as well, and it is known as Bypass transition. The term is due to 
the fact that this scenario bypasses the growth of two-dimensional waves and their 
secondary instability. Since this laminar-turbulent transition is observed even for 
values of the Reynolds number lower than the critical one, obtained by the modal 
stability theory, many shear fiows fall in the class of subcritical transitional fiows. 
The general scenario is the following. The transition does not occur simultaneously 
in the whole domain, but through nucleation and growth of organized patches of 
turbulent fiow, called turbulent spots, that eventually fill the space. The first 
observation was made by Emmons (1951) in a water table fiow. The PCf has 
been extensively studied in the past, probably because its zero mean advection 
speed allows easier tracking of the spots. Experimental investigations has shown 
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the existence of a threshold Reynolds number Re^ below which the spots keep a 
finite probability to relaminarize; among these, we remind the works by Daviaud 
et al. (1992) [Re^ = 370 ± 10), Tillmark & Alfredsson (1992) {Re^ = 360 ± 10), 
Hegseth (1996) {Rcc = 325). The most common experimental apparatus consists 
of a counter-translating belt driven by two rotating cylinders, the working fluid is 
water and a finite-amplitude disturbance is triggered by fiuid injection. Among 
the nonlinear direct numerical simulations, we report more recent works by Lagha 
& Manneville (2007), Duguet et al. (2010) {Re^ = 324 ± 1), Duguet et al. (2011) 
(Rec = 325). Usually an germ-like initial condition, (typically two counter-rotating 
vortices), is given in the physical plane. In these cases the authors showed that the 
turbulent region exhibits elongated fiow structures, called streaks, and that for PCf 
the spot shape is elliptical. When the transition is natural the streaky structure 
remains, the spots nucleate randomly in space and their shape is more irregular, or 
oblique bands are found (see Manneville, 2011). The typical distance between two 
streaks is found to be of the same order of magnitude of the channel half- height, 
Xg ^ 0{h) — 0{3h). In addition, other typical characteristics of the spot are its 
propagation speed and spreading rates which depends on the base fiow and the 
Reynolds number. However, for channel fiows the spreading of a turbulent spot 
is quite rapid, if compared to the typical turbulent diffusion: this mechanism is 
known as "growth by destabilization" of the surrounding laminar fiow. Analyzing 
the structure of the Couette spot, Lundbladth & v. Johansson (1991) and Dauchot 
& Daviaud (1995) classified it as a case between the Poiseuille and the boundary 
layer spot. 

About the latter, a boundary layer spot is characterized by a horseshoe structure. 
The complete process of transition on a fiat plate with zero pressure gradient has 
been subject to extensive studies since the beginning of the past century with 
the work of Burgers (1924) and successively by Tollmien (1929) and Schlichting 
(1933). Several features distinguish the Bbl transition from the one occurring in 
internal fiows. From the leading edge of the plate, as x increases, the laminar 
fiow is destabilized until the transition zone is reached, where arrowhead turbulent 
spots appear. A complete description of the structure and the evolution of spots 
can be found in the experimental works by Cantwell et al. (1978) and Gad-El-Hak 
et al. (1981). The former also provided beautiful visualizations of both the lateral 
side of the spot and its bottom side (the sublayer) taking advantage of the glass 
walls of the water channel. For a complete description of all the boundary layer 
mechanisms of transition and the onset of turbulence see the review by Kachanov 
(1994). 

Wave packets and role of the linear stages in the transition process 

Comparing to the amount of studies about the non-linear stages of transition 
and the descriptions of the turbulent spots, few investigations are found about the 
role of the linear evolution of small disturbances in the transitional process. The 
results of the modal analysis have probably been overestimated, and only recently a 
renewed interest in the transient evolution of linear three-dimensional disturbances 
arised. The importance of three-dimensionality and so the spanwise variation of the 
velocity components was firstly pointed out by the boundary layer experiments by 
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Klebanoff et al. (1962). The instability of this obhque wave develops in A-vortices 
(K-transition). Zang & Krist (1989) demonstrated that the growth of the oblique 
waves is correlated with the existence of a mode with a = 0, i.e. an orthogonal 
mode. The fact that the presence of this mode is a prerequisite for the rising of 
secondary instability was confirmed in earlier investigations. Schmid & Henning- 
son (1992) looked at small amplitude wave pairs, and Henningson et al. (1993) 
devoted to investigating a possible mechanism for bypass transition, pointing out 
the role of the linear phase and arguing that the mechanism for energy trans- 
fer is primarily linear. In fact, the disturbances with no streamwise dependence 
{a = 0) are usually those which experience the most rapid growth (see Fig. 4.23- 
4.26). Henningson et al. (1994) argued the necessity of linear growth mechanisms 
for subcritical growth of arbitrary amplitude perturbations. Indeed, almost all the 
Fourier components are contained in a generic initial condition, and those corre- 
sponding to the spanwise wavenumber axis are found to be rapidly excited due to 
linear mechanism. Moreover, even if the initial condition is poor in those com- 
ponents, rapid growth still occurs when non-linear interactions transfer energy in 
that area of the wavenumber space. The same happens for finite amplitude dis- 
turbances: Henningson et al. (1993) pointed out that the energy growth is only 
caused by the linear mechanism, leading to the streaky horizontal velocity pattern. 
For subcritical fiows this means that the transient growth effect must operate for 
transition to take place. 

The streaky structure, typical of various flow conflgurations, is also related to 
the evolution of optimal (in a linear sense) disturbances, which can arise and bring 
to nonlinearity (see e.g. Brandt et al., 2003). In fact, the wall-normal shape of 
linearly optimal disturbances determined by Andersson et al. (1999) is surprisingly 
similar to the measured Urms values. We should also cite the work by Cherubini 
et al. (2010), who looked for optimal initial conditions and also made a comparison 
between a linear and a nonlinear analysis of a spot evolution in Bbl. They shows 
that the streaky structure and the general shape of the spot are already determined 
by the linear analysis, due to the kinetic energy transient growth. Only in the 
following nonlinear phase, secondary instability of the streaks occurs and the spot 
central region becomes turbulent. 
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5.2 Linear spot in Plane Couette flow 

The results of the hnear superposition of a large number of waves are shown 
in Fig. 5.1-5.4. The purpose of these visualizations is to confirm the role of the 
linear transient dynamics in the complex transitional scenario, showing that in 
the evolution of a wave packet some of the typical features of a transitional flow 
may be encountered. Differently from the works by the cited authors, here a lo- 
calized disturbance is simply obtained by a superposition, with zero phase-shift, 
of a large number of waves with obliquity angle spanning the full circle. The 
polar wavenumber k is restricted to a few values chosen accordingly to the exper- 
imental evidences found in literature. For Plane Couette flow the chosen values 
are k = {5.7, 6.5, 7.3}, while the Reynolds number is 500, with reference to the 
experimental work by Hegseth (1996). Both odd and even initial condition (the 
same introduced and used in the previous chapters) are considered. The solutions 
to the initial value problem (3.1)-(3.2) are inverse-transformed, according to the 
relations (4.11) and (4.12), and then superimposed, so the complete flow fleld in an 
arbitrary domain in x and z directions can be easily obtained. If an in-phase su- 
perposition of a large number of waves is considered, a bump-like initial condition 
in the physical space is obtained, as shown in Fig. 5.2. Clearly, if the considered 
domain is wide enough, a repetitive periodic scheme can be observed. The equiv- 
alence with a two-dimensional Fourier transform is straightforward. Hence, this 
procedure is a simple way to qualitatively represent a localized perturbation in 
the physical space, containing the contributions of all obliquity angles. Moreover, 
correlations with the transient evolution of single waves, shown in Chapter 4, are 
possible. 

In the following, the results are reported in terms of the three components of 
velocity u, v, w, and the pointwise kinetic energy e{x,y,z,t) = .5 * {u'^ + w"^ + w'^) . 
Remind that the amplitude of the initial perturbation does not have any influence 
on the results, since the analysis is linear. That is the reason why all the reported 
flelds are normalized to the maximum value gained at t = 0. For Plane Couette 
flow is convenient to visualize the quantities at the channel symmetry plane y = 0, 
since the mean flow is zero. The flrst stages of the evolution of the linear spot 
seem to be characterized by a dominant spanwise rate of spreading, due to the 
faster waves with /3 — )■ 0. When the orthogonal wave (a = 0) becomes dominant, 
a streaky flow structure is found: note that a negative band of the component u 
corresponds to a positive one of v. 

In order to show the evolution of the initial perturbation in the three-dimensional 
domain, the open-source Visit tool has been used. A Matlab® script has been 
written to create the VTK flle needed by Visit, and the points structured grid 
format has been used. In Fig. 5.5, isosurfaces of streamwise velocity are shown. 
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Streamwise velocity - Plane Couette flow, Re = 500 
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Figure 5.1: Visualization of longitudinal velocity u for PCf with Re — 500. Views of xz plane 
at the channel symmetry plane at different times. The evolution of a localized perturbation is 
obtained by superposition of 220 waves with polar wavenumber k = {5.7, 6.5, 7.3}, obliquity 
angle spanning the full circle, (j> e {—90°, +90°}, with both sym and asym initial conditions. 
The values are normalized with respect to the maximum at time t = 0. 
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Wall-normal velocity - Plane Couette flow with Re = 500 
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Figure 5.2: Visualization of the wall-normal velocity v for PCf with Re = 500. Views of xz 
plane at the channel symmetry plane at different times. The evolution of a localized perturbation 
is obtained by superposition of 220 waves with polar wavenumber k = {5.7, 6.5, 7.3}, obliquity 
angle spanning the full circle, € {—90°, -1-90°}, with both sym and asym initial conditions. 
The values are normalized with respect to the maximum at time t = 0. 
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Spanwise velocity - Plane Couette flow with Re = 500 
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Figure 5.3: Visualizations of the spanwise velocity w for PCf with Re — 500. Views of xz plane 
at the channel symmetry plane at different times. The evolution of a localized perturbation is 
obtained by superposition of 220 waves with polar wavenumber k — {5.7, 6.5, 7.3}, obliquity 
angle spanning the full circle, (j) € {—90°, -1-90°}, with both sym and asym initial conditions. 
The values are normalized with respect to the maximum at time t = 0. 
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Kinetic energy - Plane Couette flow with Re = 500 
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Figure 5.4: Visualization of the kinetic energy e for PCf with Re = 500. Views of xz plane 
at the channel symmetry plane at different times. The evolution of a localized perturbation is 
obtained by superposition of 220 waves with polar wavenumber k — {5.7, 6.5, 7.3}, obliquity 
angle spanning the full circle, cj) € {—90°, +90°}, with both sym and asym initial conditions. 
The values are normalized with respect to the maximum at time t = 0. 
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3D visualization of streamwise velocity, PCf with Re = 500 
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Figure 5.5: 3D visualization of streamwise velocity for PCf with Re = 500. (a) Initial condition, 
t(j — 0; orange surface: u/uq „^ax = 0.5; blue surface: u/uq „iax — ~0.5; (b) to = 5; orange 
surface: u/uq max = 1; blue surface: u/uq max = -1; (b) to = 15; orange surface: u/uq max = 1; 
blue surface: u/uq max = ^1- Here mq max is the niaxiniuni value in the whole 3D domain at 
the initial time {uq max — 20.8). Remind that the upper wall moves in the x direction, and 
the lower wall moves in the opposite direction. As can be noticed, the initial perturbation is 
stretched, resulting in a streaky structure which is inclined with respect to the symmetry plane 
of the channel. 
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5.3 Linear spot in Blasius boundary-layer flow 

Concluding the present work, the evolution of a localized perturbation in boundary- 
layer flow is shown. The base flow here considered is the one corresponding to a 
fiat plate with zero incidence (Blasius boundary- layer, see e.g. Schlichting, 1979; 
Rosenhead, 1963). The chosen value for Reynolds number, defined with the dis- 
placement thickness is 1000, while five values of the polar wavenumber are consid- 
ered, k = {1.26, 1.57, 2.09, 3.14, 6.28}. Remind also that the spatial coordinates 
are here normalized with the displacement thiclcness, while the reference velocity 
is the free stream velocity Uoo- Simulations have been performed with two differ- 
ent initial conditions in order to get a wide database with a variety of transient 
behaviors, whose expressions are the following 

{;W=y2g-2/^ #=0 (5.1) 

-(2) 2 -v^ ■ f \ -(2) n 

Vq = y e ^ sm[Tiy) t/q = 

The former is always positive while the latter is oscillating, they both satisfy 
the boundary conditions and have their maximum near the wall, inside the bound- 
ary layer. Also for this case, the evolution of a localized disturbance obtained by 
in-phase superposition is shown. The total number of considered waves is 365. 
However, some trials have been made with a smaller number of waves, randomly 
chosen, leading to the same general conclusions. In these cases, a more irregular 
shape of the spot is observed. The affinity of the shape acquired by the wave 
packet with the one of a turbulent spots, is noticeable (see e.g. Fig. 5.6): the 
initial disturbance evolves elongating mainly in the streamwise direction, and a 
A-structure can be clearly observed. Even from noisy or dynamic initial condition 
cases, carried out by random waves superposition or random inputs in time (not 
presented in this work), it is possible to observe a flow field dominated by not 
exactly rectilinear streaks, and often a A - pattern can be recognised. 
Also in the case of Bbl, the origin of the x and z axis for non-dimensional coordi- 
nates is considered to be the location of the initial disturbance (see e.g. Fig. 5.6a). 
The three-dimensional evolution of the linear spot can be observed from Fig. 5.10, 
where isosurfaces for the streamwise velocity are shown. The qualitative behaviour 
is in agreement with the one recently shown by Cherubini et al. (2010). 
Finally, in Tab. 5.1 we report the results for a dimensional case with Uoo = 15m/s 
and V = 1.45 ■ 10^^ (air flow). This is helpful for understanding the true order of 
magnitude of the quantities involved. 
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Streamwise velocity - Blasius boundary layer flow with Re = 1000 




Figure 5.6: Visualizations of the longitudinal velocity u for Bbl with Re = 1000. Views of xz 
plane at j/o = 1-5- The evolution of a localized perturbation is obtained by superposition of 365 
waves with polar wavenuniber k — {1.26, 1.57, 2.09, 3.14, 6.28}, obliquity angle spanning the 
full circle, (p e {—90°, +90°}, and two different initial conditions. Remind that the mean flow 
moves in the right direction. 
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Wall-normal velocity - Blasius boundary layer flow with Re = 1000 




Figure 5.7: Visualizations of the wall-normal velocity v for Bbl with Re = 1000. Views of xz 
plane at yo = 1.5. The evolution of a localized perturbation is obtained by superposition of 365 
waves with polar wavenumber k — {1.26, 1.57, 2.09, 3.14, 6.28}, obliquity angle spanning the 
full circle, (p G {^90°, +90°}, and two different initial conditions. 
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Spanwise velocity - Blasius boundary layer flow with Re = 1000 




Figure 5.8: Visualizations of the spanwise velocity w for Bbl with Re — 1000. Views of xz 
plane at yo = 1-5. The evolution of a localized perturbation is obtained by superposition of 365 
waves with polar wavenumber k — {1.26, 1.57, 2.09, 3.14, 6.28}, obliquity angle spanning the 
full circle, cf) e {—90°, +90°}, and two different initial conditions. 
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Kinetic energy - Blasius boundary layer flow with Re = 1000 



30 

20 

10 

Z 

-10 

-20 

-30 

30 

20 

10 

Z 

-10 

-20 

-30 
-60 

30 

20 

10 

Z 

-10 

-20 

-30 

I 

30 

20 

10 

Z 

-10 

-20 

-30 



• ■• 



-80 



-60 -40 



-20 



20 



40 



60 



80 



^0 



-20 



20 



40 



60 



80 



100 



I 







-.* 














^ 










u 















I 

I 



1 

0.8 
0.6 
0.4 
0.2 



to 



o 




40 

30 
20 
10 



IQJ 



O 

lO 



g 

o 
»a) 



iCD 



O 



] 


20 


40 


60 


80 




100 


120 


140 






160 


18 


■ 


■ 




^^^^1 


^^^H 


^^^^^^H 


■ 




■ 




■ 


w 




1 




■ 


■ 


■ 


^^M 


^H 


^^^H 


■ 


■ 




■ 


■ 


■ 


■ 




■ 




■ 


■ 


1 




I 


■ 


■ 


1 


I 




I 


1 


H 


! 




I 




1 


1 








1 


I 


1 


1 


1 




1 


1 


1 




1 




1 


1 








1 


1 


1 


1 


■ 




■ 


■ 


H 


1 




■ 




■ 


■ 








I 


■ 


1 
















■ 








■ 




m 


^^^^^^1 


^^1 


^^^^^^H 






^^1 




^H 










■ 


^^1 


■ 



25 
20 

15 
10 
5 



o 



15 



(03 



O 
iO 



10 



60 



80 



100 



120 



140 

X 



160 



180 



200 



220 



I 



o 



)QJ 



Figure 5.9: Kinetic energy visualizations for Bbl with Re = 1000. Views of xz plane at yo — 
1.5. The evolution of a localized perturbation is obtained by superposition of 365 waves with 
polar wavenuniber k = {1.26, 1.57, 2.09, 3.14, 6.28}, obliquity angle spanning the full circle, 
(p E {—90°, +90°}, and two different initial conditions. 
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5.3 - Linear spot in Blasius boundary-layer Row 



3D visualization of streamwise velocity, Bbl with Re = 1000 




180 



Figure 5.10: 3D visualization of streamwise velocity for Bbl with Re = 1000. (a) Initial con- 
dition, to — 0; orange surface: u/uq „iax — 0.5; blue surface: u/uq max =^ —0.5; (b) ip = 50; 
orange surface: u/uq max = 1; blue surface: u/uq max = ^1; (b) to — 150; orange surface: 
u/uo max = 1; blue surfacc: u/uq max = — 1- Here uq max is the niaxiniuni value in the whole 3D 
domain at the initial time (uq max — 214). 
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5 - Wave packets linear evolution 



In the following table, a summary of the linear spot evolution for the analyzed 
case is shown, and also the dimensional quantities are reported for a specific case 
with Uoo = Ibm/s and air flow (z/ = 1.45 ■ 10~^). In this conditions, a Reynolds 
number of 1000 is found at a distance x — xie = 328 mm from the leading edge 
of the fiat plate, where 6* = 0.97 mm. This is considered the origin of the non- 
dimensional X axis, where the linear spot is triggered. Actually, in the evolution of 
the perturbation the boundary-layer thickness and the Reynolds number change, 
but in this analysis the approximation of nearly-parallel flow applies. In Tab. 5.1 
Xc is the location of the center of the spot, Uc is the longitudinal velocity of the spot 
center, L^ and L^ are the longitudinal and spanwise sizes for yo = 1.5, respectively. 
All these quantities are esteemed from visualizations of the streamwise perturbation 
velocity. 

t = t = 5 t = 50 t = 150 t = 250 
ms 0.323 ms 32.3 ms 9.70 ms 16.2 ms 



6* 


0.967 mm 


0.972 mm 


1.01 mm 


1.08 mm 


1.16 mm 


^0.99 


2.770 mm 


2.773 mm 


2.88 mm 


3.08 mm 


3.31 mm 


"^ c 



328 mm 


2 
330 mm 


29.6 
357 mm 


88.0 
413 mm 


148 
471 mm 


C/c 




m/s 


0.8 

12 m/s 


0.67 

10 m/s 


0.64 

9.6 m/s 


0.64 

9.6 m/s 


L^ 


6 

5.82 mm 


12 
11.6 mm 


25 
24.3 mm 


60 

58.2 mm 


94 
91.2 mm 


L. 


6 

5.82 mm 


12 
11.6 mm 


18 
17.5 mm 


24 
23.3 mm 


24 
23.3 mm 



Table 5.1: Dimensional quantities for the evolution of a localized perturbation in a Blasius 
boundary-layer flow with Re = 1000. The approximation of nearly-parallel flow is applied, in 
fact the Reynolds number for every simulation is fixed. It can be noticed that the true spatial 
and temporal scales are quite small. 
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Chapter 6 



Conclusions 



The present work deals with the hydrodynamic stabihty non-modal analysis. 
In the first part, a close form solution to the three-dimensional Orr-Sommerfeld 
and Squire IVP in the form of orthogonal functions expansion was researched. The 
Galerkin variational method was then successfully implemented in the Matlab® en- 
vironment to numerically compute approximate solutions to the coupled equations, 
for bounded flows. The Chandrasekhar functions revealed to ensure a convergence 
rate scaling as A^^ to the correct solution even for the non-modal analysis. The 
advantages of this method can be summarized in the independence of the accuracy 
on both the temporal and spatial grids, which can be considered arbitrary, and the 
very low time computational cost. Moreover, since there is no "marching" in time, 
no stiffness problems are encountered and accurate solutions can be obtained up 
to very high times. The spectra of the Orr-Sommerfeld and Squire operators are 
computed with high precision, as well. 

The code has been intensively used to focus on the temporal evolution of the 
wave frequency and phase velocity, poorly investigated in the past. The results 
confirmed recent observations about the frequency jump in the v component of 
flow velocity, considered as the end of the Early transient. After this first jump, 
the frequency of v for Plane Couette flow experiences a periodic modulation about 
the asymptotic value, which has been motivated and investigated in detail. 

A new result is the presence of a second jump in the phase velocity of the vorticity 
component fj and consequently of the other components of velocity, typically for 
high times. The presence of a second jump and the possibility for different values 
of asymptotic frequency of the signals were motivated and no contradictions with 
the results of the modal theory subsist. This is the proof of the existence of an 
Intermediate transient., in fact, only after the last jump the solution reaches its 
asymptotic state. Moreover, a connection between the frequency jumps and the 
establishing of a self-similarity condition in time for both the velocity and vorticity 
profiles was found and investigated for both Plane Couette flow and Plane Poiseuille 
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6 - Conclusions 



flow. The behavior in the physical space was also shown. 

The last result deals with the linear evolution of wave packets. Through super- 
position of waves with limited wavenumber range a wave packet is reconstructed 
for Plane Couette flow and Blasius boundary-layer flow. The structure of the linear 
spot revealed to have many common features with the early stages of a turbulent 
spot, particularly the streaky structure and the shape. This is in agreement with 
recent ideas and observations and supports the thesis of the underrated importance 
of the linear mechanisms such as the transient growth, in the transitional scenario. 
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Appendix A 



A.l The basis eigenfunctions 

The solutions of the problem (3.9) with homogeneous boundary conditions (3.10) 
can be found by seeking a solution of the following form 

X{y) = Cie^y + C2e-^y + C^e'^y + de-'^y (A.l) 

where Cj are constants. Since we are interested in real eigenfunctions, the above 
expression becomes 

X{y) = Asinh{Xy) + Bcosh{Xy) + Csin{Xy) + Dcos{Xy) (A.2) 

where A, B, C and D are constants, to be determined. Substituting the solution 
in the expressions of the boundary conditions, a system of four algebraic equations 
where the constant are the unknown is found. The following relations are obtained 

stnh[X) cosh{X) 



(A.4) 



- tanh{x) + c^'^(^) -cos{X)tanh{X) - sin{X) l'C\ _ fO 
-iSS) + "^^(^) cos{X)tanh{X) + stn{X) \ [dJ " ^0 

Nontrivial solution to the system (A.4) are possible if the determinant of the coef- 
ficients matrix vanishes, this leads to the equation 

tan'^{X)-tanh^{X) = (A.5) 

from which the couple of trascendental equations (3.11) and (3.12) are derived. The 
solution to these equations consists in two sets of infinite and discrete eigenvalues 
A„, that can be computed numerically by the bisection or the Newton- Raphson 
method (Fig. A.l). In Tab. A.l the first twenty eigenvalues are shown. 
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Figure A.l: Graphic solution of the eigenvalue problem (A. 4). Odd indices indicate the eigen- 
values corresponding to odd eigenfunctions and even indices correspond to even eigenfunctions. 
Since both sets are antisymmetric with resprect to A = axis, only positive eigenvalues are 
considered. 



n 


An {odd) 


n 


'^n (even) 


1 


3.926602 


2 


2.365020 


3 


7.068582 


4 


5.497803 


5 


10.21017 


6 


8.639379 


7 


13.35176 


8 


11.78097 


9 


16.49336 


10 


14.92256 


11 


19.63495 


12 


18.06415 


13 


22.77654 


14 


21.20575 


15 


25.91813 


16 


24.34734 


17 


29.05973 


18 


27.48893 


19 


32.20132 


20 


30.63052 



(a) Odd eigenvalues 



(b) Even eigenvalues 



Table A.l: First 20 eigenvalues, numerically computed using the Newton-Raphson method. 
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Appendix B 



B.l Matlab® scripts for channel flows Galerkin method 
B.1.1 Main program: "main_ ivp_ galerkin. m" 

1 %% IVP SOLUTION BY GALERKIN METHOD for Channel flows 

2 % The stability IVP is reduced to an ODE system and solved by 

3 % Galerkin variational method. 

4 clear aU , close all , clc 

5 tic 

6 %% Simulation parameters 

7 global tol N U DU D2U type Re k_polar alfa beta phi ss y y_ic t 

8 type=' Couette ' ; 

9 Re = 500; % Reynolds number 

10 k_polar=6.5; % Polar wavenumber 

11 phi=5; % Obliquity p erturb ation angle 

12 alf a=k_polar*cosd ( phi ) ; 

13 beta=k_polar*sind ( phi ) ; 

14 ss='syni'; % Type of initial condition 

15 h = 0.01; 

16 y= — l:h:l; % Grid for solution 

17 y_ic = — l:le — 4: 1 ; % Fine grid for initial condition computaiton 

18 N=200; % Num. of eigenmodes chosen for calculation (must be even). 

19 tol=le — 10; % Tolerance for error checking . 

20 col=colormap( jet ( 1 5) ) ; 

21 y0=0.5; 

22 temp=find(y>=yO); 

23 idxy=temp ( 1 ) ; 

24 dt=0.5; 

25 t=0:dt:100; % Time grid 
26 

27 disp 
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28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 



48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
70 
71 
72 
73 
74 
75 
76 
77 
78 
79 
80 
81 
82 
83 



disp ( [ ' S i mulat ion ^ of „ Plane „ ' ,nuin2str(type) , '„flow ' ] ) 

disp ( [ 'Re=^ ' ,num2str (Re) , ' , „k_polar=„ ' ,num2str (k_polar ) , ' , „phi=„ 

nuin2str ( phi ) , '^deg , ^ ' ,nuin2str ( ss), '^Initial^Condition. ']) 
disp ( [ 'N=^ ' ,num2str (N) , '^eigen functions ']) 

disp mmmmmmmmmmmmmmmmmmmmmmmmmmmmwm: 



%% Base flow definition 

if strcmp(type, ' Couette ')==1 

U=y; 

DU=1; 

D2U=0; 
elseif strcmp(type, 'Poiseuille ')==1 

U=(l-y.-2); 

DU=-2*y ; 

D2U=-2; 
end 

% 

% 

%% Eigenmodes computation 

% The chosen base of functions is composed by two set of eigenfuncions of 
% the problem X_yyyy=lam'^4*X with b.c. X(- 1)=X(1) =X'(- 1) =X'( 1 ) = 0. 
% One set is made of odd modes, the other contains even modes. Both sets are 
% necessary to obtain a correct solution of the 3D problem , independently 
% on the parity of the initial condition . 

% Here we put togheter the two set in following order:[lo le 2o 2e So 3e 
%(e=even, o=odd), and the corresponding eigenvalues are ordered in this 
%sequence too. See " compute_ eigenvalues .m" for their computation . 



45 
46 


% 

%% NORMAI^VELOCITY SOLUTION 


v(y 


t) 


47 


% 







1 



load ' eigenvalues . mat 
g=gtot(l:N); 



;% " eigenvalues .mat " , contains the eigenvalues 
%corresp onding to form functions . 
% g are the eigenvalues. 



% first derivatives of eigenmodes 
% second derivatives of eigenmodes 



eigenmodes=zeros ( length (g) , length (y ) ) ; 
Dleigenniodes=zeros ( length (g) , length (y ) ) 
D2eigenniodes=zeros ( length (g) , length (y ) ) 
for i =l:length(g) 

if niod(i,2)==l % odd eigenfuncion 
eigenmodes (i ,:) = l/sqrt (2) * ( sinh(g( i )*y)/ sinh (g( i )) -sin (g( i )*y)/ sin (g ( i ))); 
eigenmodes_ic(i ,:) = l/sqrt(2)*( sinh (g(i )*y_ic)/sinh(g(i)) — sin(g(i )*y_ic)/sin(g(i ) )) 



Dl eigenmodes (i ,:) = g(i)/sqrt(2)*( cosh(g ( i )*y)/sinh(g( i)) — cos(g 
D2eigenmodes (i,:) = g(i)''2/sqrt(2)*(sinh(g(i)*y)/sinh(g(i)) + sin 
elseif mod(i,2) ==0 % even eigenfuncion 
eigenmodes(i ,:) = l/sqrt(2)*( cosh(g( i )*y)/cosh(g( i )) — cos (g( i )*y 



eigenmodes_ic(i ,:) = l/sqrt(2)*( cosh(g( i )*y_ic)/cosh(g( i)) — cos(g(i)*y_ic)/ cos (g( i ) ) ) 



Dl eigenmodes (i ,:) = g(i)/sqrt(2)*( sinh (g ( i ) *y )/cosh(g( i )) + sin (g 
D2eigenmodes (i ,:) = g(i)''2/sqrt(2)* (cosh(g( i ) *y)/cosh(g( i )) + cos 
end 
end 



i)*y)/sin(g(i ))); 
g(i)*y)/sin(g(i ))) 

/cos(g(i ))); 



i)*y)/cos(g(i ))); 
g(i)*y)/cos(g(i ))) 



%% Computation of c o efficients CnO , from pres crib ed initial condition 
if strcnip(ss , 'sym')==l 

ic=(l-y.-2).-2; 

ic_a=(l— y_ic . " 2 ) . " 2 ; % symmetric initial condition 

elseif strcnip(ss , 'asym')==l 
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B.l - Matlab® scripts for channel Hows Galerkin method 



84 ic=(l-y.-2).-2; 

85 ic_a=y .*(1 — y_ic . " 2 ) . " 2 ; % asymmetric initial condition 

86 end 
87 

88 A=zeros(l ,N); 

89 approx_ic=zeros ( 1 , length (y )) ; 

90 CnO=zeros(l ,N); 

91 for i=l:N; 

92 num=trapz (y_ic , ic_a . * eigenniodes_ic ( i ,'■))', 

93 den = l; 

94 CnO( i)=nuin/den ; 

95 approx_ic (1 ,:)= approx_ic (1 , :) +CnO( i )* eigenmodes ( i , : ) ; 

96 end 
97 

98 %% Computation of matrices D S F Ul U2 US H G 

99 % ODE system is [H] { dc/dt}-[G] { c}=0 

100 D=zeros(N,N); 

101 S=zeros(N,N); 

102 I^zeros(N,N); 

103 F=zeros(N,N); 

104 Ul=zeros(N,N); 

105 U2=zeros(N,N); 

106 G=zeros(N,N); 

107 H=zeros(N,N); 
108 

109 %% D and F 

110 Dnm=zeros (1 ,N) ; 

111 Fnin=zeros (1 ,N) ; 

112 for i=l:N 

113 nDm(i) = l; 

114 Fnm( i)=Dnm( i )*g( i ) "4; 

115 end 

116 D=diag(Dnm,0); 

117 F=diag(Fnm,0); 
118 

119 bbb=(cosh(2*g)-cos(2*g)) . / ( sinh (2*g)-sin (2*g ) ) ; 

120 idNaN=find(isnan(bbb) = = l); 

121 bbb(idNaN) = l; 

122 S=zeros(N,N); 
123 

124 for ni=l:N 

125 for n=l:N 

126 %% S 

127 if mod((m^n),2) = = &fe n~^n 

128 S(m,n)=+4*g(n)'-2*g(m)'-2/(g(n)~4-g(ni)~4)*(g(n)*bbb(n)-g(ni)*bbb(m)) ; 

129 elseif mod( (irH-n) ,2) = = 1 &fe n~^n 

130 S(m,n) = 0; 

131 elseif n=^ii 

132 S(m,n) = -(g(n)^2*bbb(n)~2-g(n)*bbb(n)); 

133 end 
134 

135 %% Ul and U2 

136 if strcnip(type , ' Couette ')==1 % analytical expressions for Gouette 
137 

138 if mod((m^n),2) = = l &fe n~^n 

139 Ul(m,n)=4*g(n)'-2*g(m)'-2/(g(n)'-4-g(ni)~4)*(-l + g(n)*bbb(n) . . . 
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140 -g(m)*bbb(m))-8*(g(n)"4+g(m)"4)/(g(n)"4-g(m)"4)"2*g(n)"2*g(m)"2; 

141 end 

142 if niod((m^n),2) = = l &fe n~^n 

143 U2(m,n) = 16*g(n)~3*g(m)'^3*bbb(n)*bbb(m)/(g(n)'^4-g(m)'^4)~2; 

144 elseif mod( (nifn) ,2) = = || n=^ii 

145 U2(m,n)=0; 

146 end 

147 elseif strcnip(type , ' Poiseuille ')==1 % numerical integration 
148 

149 prod=U. * D2eigenniodes (n , : ) . * eigenmodes (m, : ) ; 

150 Ul(m, n)=trapz (y , prod ) ; 

151 clear prod 

152 prod=U. * eigenmodes (n , : ) . * eigenmodes (m, : ) ; 

153 U2(m, n)=trapz (y , prod ) ; 

154 end 
155 

156 end 

157 end 

158 %% US 

159 U3=D2U*D ; 
160 

161 %% H and G 
162 

163 H=S-k_polar^2*D; 

164 G=-li*alfa*Ul. .. 

165 +li * alf a*k_polar "2*U2 . . . 

166 +li*alfa*U3. . . 

167 +l/Re*F . . . 

168 -2*k_polar '^2/Re*S . . . 

169 +k_polar'^4/Re*D; 
170 

171 %% ODE system solution 

172 % Compute the series. The solution V has length (y) rows and length (t) 

173 % colum,ns . 

174 % V(y,t)=sum( Cn(t)*Xn(y)) 
175 

176 A=H\G; 

177 % Check error in H matrix inversion % 

178 relerrA_norm2=norm( (H*A— G) ,2) /norm(G, 2 ) ; 

179 relerrA_norminf=norm( (H*A— G) . inf )/norm(G, inf ) ; 

180 disp ( ' ' ) 

181 disp ( [ ' Absolute^error „of „A„computation„using„\„in„norm„ 2: „ ' ,... 

182 num2str (relerrA_norm2*norm(G, 2 ) ) ] ) 

183 disp ( [ ' Relative^error ^of „A„computation^using„\^in„norm„ 2: „ ' ,... 

184 num2str (relerrA_norm2 ) ] ) 

185 disp ( [ ' Relative^error ^of ^A^computation„using„\„in„norm„Inf : „ ' ,... 

186 num2str (relerrA_norminf ) ] ) 

187 disp(' ') 

188 % % 

189 

190 % Diagonalizing matrix A 

191 [ L, LAMBDA] = eig (A); 
192 

193 Cn=zeros(N,length(t )) ; 

194 hh=zeros (N, length(t ) ) ; 

195 hO=L\(CnO'); 
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B.l - Matlab® scripts for channel Hows Galerkin method 



196 

197 % Check error in L matrix inversion 

198 relerrL_norni2=norm( (L*hO— CnO ' ) , 2) /norin(CnO , 2 ) ; 

199 relerrL_norniinf=norm( (L*hO— CnO ' ) , inf )/norm(CnO , inf ) ; 



200 disp ( ' ' ) 

201 disp ( [ ' Absolute^error ^of „hO„coniputation„using„\„in„norni„ 2: „ ' ,... 

202 nuni2str (relerrL_norm2*norni(CnO , 2 ) ) ] ) 

203 disp([ ' Relative^error „of „hO„computation„using„\„in„norni„ 2 : „ ' ,... 

204 num2str (relerrL_norni2 ) ] ) 

205 disp([ ' Relative ^ error „of„hO„coinputation„using„\„in„norni„ Inf :„ ' ,... 

206 num2str ( relerrL_norniinf ) ] ) 
207 

208 if relerrL_norni2>tol 
209 

210 [h02 ,flag_hO ,relerrL2_norni2] = gmres (L , CnO ' , [] , tol ,N) ; 

211 if ( relerrL2_norm2) >(relerrL_norni2 ) 

212 disp ([ 'Won' 't„useXME«. ']) 

213 elseif flag_hO==0 &fe (relerrL2_norni2) <(relerrL_norni2 ) 

214 disp ( [ 'GMRS^converged^ within „ tolerance • ^ ' ] ) 

215 disp ( [ ' Absolut e„ err or „ of „hO„computation „ using „GMFiS„ in „norni„ 2 : „ ' , . . . 

216 num2str(norm((L*h02-CnO') ,2))]) 

217 disp ( [ ' Relative^error ^of ^h0^coniputation„using„GMFiS„in„norni„ 2 : „ ' , . . . 

218 nuin2str ( relerrL2_norm2 ) ] ) 

219 disp([ 'Relative^error „of „liO„coniputation„using„\„in„norni„ Inf : „ ' , . . . 

220 num2str ( relerrL_norniinf ) ] ) 

221 h0=li02 ; 

222 end 

223 end 

224 disp(' ') 

225 

226 

227 % Compute coefficients Cn(t) 

228 lili=zeros(N,length(t )) ; 

229 for i=l:N 

230 lili(i ,:)=liO(i)*exp(LAMBDA(i , i)*t ); 

231 end 
232 

233 Cn=L*(hli); % Coefficients of the series Cn(t) (N rows , length ( t ) columns) 
234 

235 % Compute v(y,t), v_y(y,t), v_yy(y,t) 

236 v=zeros ( length (y ) , length ( t ) ) ; 

237 Dv^zeros (length (y ) , length ( t ) ) ; 

238 D2v=zeros(length(y) ,length(t )) ; 

239 for j=l:length(y) 

240 for i=l:N 

241 v(j ,:) = v(j ,:)+ Cn( i , : ) * eigenmodes ( i , j ) ; 

242 Dv(j ,:)=Dv(j ,:) + Cn(i , : ) * Dleigenmodes ( i , j ); 

243 D2v(j ,:)=D2v(j ,:)+Cn(i , : ) * D2eigenmodes ( i , j ); 

244 end 

245 end 
246 

247 % Check error in initial condition computation % 

248 abserr_v_norni2=norni(v(: , 1) — ic ' ,2); 

249 abserr_v_norniinf=norni(v(: ,1)— ic',inf); 

250 disp(' ') 

251 disp([ 'Absolute^ error „on„v( t=0)^in^norni^ 2: ^ ' ,nuni2str (abserr_v_norni2 ) ] ) 
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252 disp([ 'Absolute^ error „ on „v( t =0)„in„norni„ inf : „ ' ,nuin2str ( abscrr_v_norniinf ) ] ) 

253 disp ( ' ' ) 

254 

255 % % 

256 %% NORMAI^VORTICITY SOLUTION 

257 % % 

258 % calls function solve _ squire .ni 

259 % oy(y,t) is omega_y (eta) 

260 oy=solve_ squire (L, LAMBDA, hO , eigenmodes , g ) ; 
261 

262 % % 

263 %% POSTPROCESSING 

264 % V and eta modulus at all y and times 

265 niod_v=abs ( v ) ; 

266 niod_oy=abs ( oy ) ; 

267 %% Streamwise velocity u(y,t) 

268 u = l/(li*k_polar" 2)*(beta*oy— alf a*Dv) ; 

269 niod_u=abs (u ) ; % modulus at all y and all times 

270 wrphase_u=phase (u( idxy , : ) ) ; 

271 phase_u=unwrap( wrphase_u ) ; 
272 

273 %% Spanwise velocity w(y,t) 

274 w=l/(l i *k_ polar "2)*( — beta*Dv— alfa *oy ) ; 

275 mod_-w^abs (w) ; 

276 wrphase_w=phase (w( idxy , : ) ) ; 

277 phase_w=unwrap(wrphase_w ) ; 

278 %% Frequency and phase velocity 

279 % Calculated from v (4th order finite differences centered scheme) 

280 wrphase_v=phase (v(idxy , : ) ) ; 

281 phase_v=unwrap( wrphase_v ) ; 

282 jj =3:length(t)-2; 

283 omega_v=zeros ( 1 , length( t ) ) ; 

284 omega_v( j j )=(+phase_v ( jj —2)— 8* phase _v ( jj --l)+8*phase_v( j j + 1) . . . 

285 -phase_v(jj +2))/(12*dt ); 

286 omega_v(l) = ( — 25*phase_v(l)+48*phase_v(2) — 36*phase_v(3) + 16*phase_v ( 4) . . . 

287 -3*phase_v(5))/(12*dt); 

288 omega_v(2) = (— 3*phase_v(l) — 10*phase_v(2) + 18*phase_v(3) —6* phase _v (4) + . . . 

289 phase_v(5))/(12*dt); 

290 oniega_v(end— 1)=— (— 3*phase_v (end) — 10*phase_v(end— l) + 18*phase_v (end — 2) . . . 

291 -6*phase_v(end-3)+phase_v(end-4))/(12*dt ) ; 

292 omega_v(end)=—(— 25* phase _v (end)+48*phase_v(end — 1)— 36*phase_v (end — 2) . . . 

293 +16*phase_v(end-3)-3*phase_v(end-4))/(12*dt ) ; 
294 

295 c_v=abs (omega_v)/ k_polar ; % Phase velocity 

296 

297 

298 %% Energy Crowth Factor C = integral (u2+v2+w2) in [-1 1] 

299 ke=niod_u. ~2+niod_v. ~2+mod_w. " 2 ; % kinetic energy in wavespace 

300 eO=trapz(y,ke (: ,1)); 

301 G=trapz(y,ke,l)/eO; 
302 

303 %% PLOTS 

304 toe 
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B.1.2 Function: "solve_ squire. m" 

1 function [ E]= squire _particular_solution(L ,LA]VIBDA, hO , eigenmodes , g) 
2 

3 global tol type U DU N Re k_polar alfa beta y t 

4 %% Basic eig enf unctions computation 

5 j=l:N/2; 

6 g_sq_odd=j *pi ; 

7 g_sq_even = (2*j — l)*pi /2; 



9 g_sq=zeros (1 ,N) ; 

10 g_sq(l:2:N)=g_sq_odd; 

11 g_sq (2 : 2 :N)=g_sq_even ; 
12 

13 eigenmodes _squire=zer OS (length (g_sq) , length (y ) ) ; 

14 D2eigenmodes_squire=zeros (length (g_sq) ,length(y)); 

15 for i =1 :2: length(g_sq) 

16 eigenmodes _ squire (i ,:) = sin(g_sq(i)*y); 

17 D 2 eigenmodes _ squire (i ,:) = — g_sq(i)"2*sin (g_sq( i )*y ) ; 

18 end 

19 for i=2:2:length(g_sq) 

20 eigenmodes_squire(i,:) = cos (g_sq(i)*y); 

21 D 2 eigenmodes _ squire (i,:) = — g_sq(i)~2*cos(g_sq(i)*y); 

22 end 
23 

24 %% Compute matrices 
25 

26 DD=zeros(N,N) 

27 SS=zeros(N,N) 

28 lJU=zeros(N,N) 

29 GG=zeros(N,N) 

30 BB=zeros(N,N): 

31 FOR^zeros(N,N); 
32 

33 %% D 

34 DD=eye(N); 

35 %% SS 

36 SS=zeros(N,N); 

37 for n = l:N 

38 SS(n,n)=-g_sq(n)^2 

39 end 

40 %% UU and FORZ 

41 for m=l:N 

42 for n=l:N 
43 
44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 



if strcmp(type , ' Couette ')==1 % Analytical expression 
if mod( (n-kn) ,2) = = 1 

lJU(m, n) = (-l)~((m^n + l)/2)*4*g_sq(m)*g_sq(n ) . . . 
/(g_sq(n)-2-g_sq(m)-2)-2; 
end 



if mod(m,2) = = l &fe mod(n,2)==l 

FORZ(m,n) = sqrt(2)*2*g_sq(m)*g(n)"2*(-l)"((m+l)/2)/(g_sq(m)"4-g(n)~4); 
elseif mod(n,2)==0 &fe mod(m,2) ==0 

FORZ(m,n) = sqrt (2)*2*g_sq(m)*g(n)~2*(-l)~(m/2)/(g_sq(m)~4-g(n)~4); 
end 
elseif strcmp(type, 'Poiseuille ')==1 % Numerical integration 
prod=U. * eigenmodes _ squire (n , : ) . * eigenmodes _ squire (m, : ) ; 
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55 UU(m, n)=trapz (y , prod) ; 

56 clear prod 

57 prod=DU. * eigenmodes (n , : ) . * eigeninodes_squire (m, : ) ; 

58 F0RZ2(ni,n)=trapz(y,prod); 

59 end 

60 end 

61 end 
62 

63 %% Find particular solution E_p(y,t) 
64 

65 GG=-li*alfa*lJU+l/Re*SS-k_polar^2/Re*DD; 

66 BB=-li*beta*FORZ*L; 
67 

68 % Compute coefficient matrix aa 

69 aa=zeros(N,N) ; 

70 idx_ginrs = 0; 
71 

72 for i=l:N 

73 TT=iAlvIBDA( i , i ) * eye (N) -GG; 

74 rhs=BB(: , i)*hO(i ); 

75 aa(: , i)=TT\rhs ; 
76 

77 % Check error in TT matrix inversion % 

78 relerrTT_norni2 ( i )=norm( (TT*aa (: , i) — rhs) ,2) /norni(rhs , 2 ) ; 

79 relerrTT_norminf( i)=norm( (TT*aa ( : , i ) — rhs ) , inf )/norni(rhs , inf ) ; 

80 abserrTT_norin2 ( i)= relerrTT_norm2 ( i )*norm(rhs , 2) ; 

81 if relerrTT_norni2(i)>tol 

82 [aa2(: , i ) , flag_TT , relerrTT2_norni2 ( i )] = gmres (TT, rhs , [] , tol ,N) ; 

83 if (relerrTT2_norm2(i))<(relerrTT_norni2(i )) 

84 aa ( : , i ) = aa2 ( : , i ) ; 

85 relerrTT_norin2 ( i )=relerrTT2_norni2 ( i ) 

86 relerrTT_norniinf(i )=norni( (TT*aa (:,i) — rhs),inf )/norni( rhs , inf ) ; 

87 abserrTT_norni2 ( i )= relerrTT_norni2 ( i )*norm( rhs , 2 ) ; 

88 idx_gnirs = [idx_gnirs i ] ; 

89 end 

90 end 

91 % % 

92 end 

93 if idx_gnirs ~=0 

94 disp ( [ ' Used „GMFiS„ for „aa , „ at „columns„ of „ indices : „ ' ,nuni2str (idx_svd ) ] ) 

95 end 

96 disp ( [ 'Max. ^absolute^error „for „aa„in„norni„ 2: „ ' , nuni2str (max(abserrTT_norni2 ) ) ] ) 

97 disp ( [ 'Max. ^ relative „ error ^for ^aa^in„norm„ 2: „ ' , nuni2str (max(relerrTT_norni2 ) ) ] ) 

98 disp ( [ 'Max. ^ relative^ error ^for ^aa^in^norm^Inf : „ ' , num2str (max(relerrTT_norminf ) ) ] ) 
99 

100 % Compute c o effi cients for particular solution 

101 expon=zeros (N, length ( t )) ; 

102 for i=l:N 

103 expon(i ,:)=exp(LATvIBDA(i , i)*t); 

104 end 

105 Jp=aa*expon; 

106 JpO=KK*ones (N, 1 ) ; % for zero initial vorticity 

107 % Compute solution 

108 E_p=zeros ( length (y) , length (t ) ) ; 

109 for j=l:length(y) 

110 for i=l:N 



94 



B.l - Matlab® scripts for channel Hows Galerkin method 



111 E_p(j ,:)=E_p(j ,:)+ Jp ( i , : ) * eigenmodes_squire ( i , j ) ; 

112 end 

113 end 
114 

115 %% Compute homogeneous solution E_h(y,t) 
116 

117 [L_eta,LAMBDA_eta] = eig(GG); 

118 hhO_eta=— L_eta\ JpO ; % Here we use initial condition eta ( t =0 ,y)=0 
119 

120 % Check error in L_eta matrix inversion 

121 relerrLeta_norni2=norni( (L_eta*hhO_eta+JpO) ,2) /norni(— JpO , 2 ) ; 

122 relerrLeta_norininf=norni( (L_eta*hhO_eta+JpO ) , inf )/norni( — JpO , inf ) ; 

123 abserrLeta_norni2= relerrLeta_norni2*norni(— JpO , 2 ) ; 

124 disp ( ' ' ) 

125 disp ( [ ' Absolute^error ^of „hh0_eta„coniputation„using„\„in„norni„2 : ^ ' ,... 

126 nuni2str ( relerrLeta_norm2*norni(— JpO , 2 ) ) ] ) 

127 disp([ ' Relative^error „of ^hh0_eta„computation„using„\„in„norni„2 : „ ' ,... 

128 nuni2str ( relerrLeta_norni2 ) ] ) 

129 disp([ ' Rclative^error ^of ^hhO_eta„coniputation„using„\„in„norni„Inf : „ ' ,... 

130 nuni2str ( abserrLeta_norm2 ) ] ) 
131 

132 if relerrLeta_norni2<tol 

133 [ hh0_eta2 , flag_hhO_eta , relerrLeta2_norni2 ] = gmres (L_eta,— JpO , [] ,tol,N); 

134 if ( relerrLeta_norm2) >(relerrLeta2_norni2 ) 

135 disp ([ 'Won' 't^use^GMRS. ']) 

136 elseif flag_hhO_eta==0 && ( relerrLeta2_norni2) <(relerrLeta_norm2 ) 

137 disp ( [ 'GMlRS„converged„ within „ tolerance . „ ' ] ) 

138 disp ( [ ' Absolute^error „of ^hhO_eta„computation„using„G]VIRS„in„norm„2: „ ' ,. 

139 nuni2str ( relerrLeta_norm2*norni(— JpO , 2 ) ) ] ) 

140 disp ( [ ' Relative^error „of „hhO_eta„computation„using^GMFiS„in„norni„2 : „ ' ,. 

141 nuni2str ( relerrLeta_norm2 ) ] ) 

142 disp ( [ ' Relative^error „of „hhO_eta„computation„using„\„in„norm„Inf : „ ' ,.. 

143 nuni2str ( relerrLeta_norniinf)] ) 

144 hh0_eta=hh0_eta2; 

145 end 

146 end 

147 % % 

148 

149 hh_eta=zeros (N, length ( t )) ; 

150 for i=l:N 

151 hh_eta(i , : ) = hhO_eta( i ) *exp(LAMBDA_eta( i , i)*t ); 

152 end 

153 % Compute coefficients for homogeneous solution 

154 Jh=L_eta*(hh_eta) ; % Size of Jh is ( N x length (t ) ) 

155 % Compute homogeneous solution 

156 E_h=zeros ( length (y) , length ( t ) ) ; 

157 for j =l:length (y) 

158 for i=l:N 

159 E_h( j , :) =E_h( j ,:)+ Jh(i ,:)*eigenmodes_squire(i ,j ); 

160 end 

161 end 
162 

163 %% Compute complete solution E(y,t) 

164 E=E_p+E_h; 
165 

166 return 
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